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ABSTRACT 

Magnetic protection of potentially habitable planets plays a central role in determining their actual 
habitability and/or the chances of detecting atmospheric biosignatures. We develop here a thermal 
evolution model of potentially habitable Earth-like planets and super-Earths. Using up-to-date dy- 
namo scaling laws we predict the properties of core dynamo magnetic fields and study the influence 
of thermal evolution on their properties. The level of magnetic protection of tidally locked and un- 
locked planets is estimated by combining simplified models of the planetary magnetosphere and a 
phenomenological description of the stellar wind. Thermal evolution introduces a strong dependence 
of magnetic protection on planetary mass and rotation rate. Tidally locked terrestrial planets with 
an Earth-like composition would have early dayside magnetospause distances between 1.5 and 4.0 i?p, 
larger than previously estimated. Unlocked planets with periods of rotation ^ 1 day are protected by 
magnetospheres extending between 3 and 8 Rp. Our results are robust against variations in planetary 
bulk composition and uncertainties in other critical model parameters. For illustration purposes the 
thermal evolution and magnetic protection of the potentially habitable super-Earths GL 581d, GJ 
667Cc and HD 40307g were also studied. Assuming an Earth-like composition we found that the 
dynamos of these planets are already extinct or close to being shut down. While GL 581d is the 
best protected, the protection of HD 40307g cannot be reliably estimated. GJ 667Cc, even under 
optimistic conditions, seems to be severely exposed to the stellar wind and, under the conditions of 
our model, has probably suffered massive atmospheric losses. 

Subject headings: Planetary systems - Planets and satellites: magnetic fields, physical evolution - 
Planet-star interactions 



1. INTRODUCTION 

The discovery of extrasolar habitable planets is one of 
the most ambitious challenges in exoplanetary research. 
At the time of writing, there are almost 861 confirmed 
exoplanetf0 including 61 classified as Earth-like planets 
(EPs, M ^ 1 Mm) and sup er-Earths (SEs, M - 1 - 
10 M®, [Valencia et "all 120061 hereafter VAL06). Among 
these low mass plane ts ther e are three confirmed SEs, GJ 

2007: 
I20T2I1. 



667Cc (Bonfils et aT1 [20TTI). GL 581 d (Udrv ct „ 
iMavoret al.ll2009| ). and HD 40307g (iTuomi et al 



and tens of Keple r candidates (jBorucki et al 



, I20T1I: 

iBatalha et al.|[2012[ ) that are close or insid e the habit- 
able z one (HZ) of thei r host st ars (see e. g. 'Selsi s et al.l 
[2007t iPepe eFaLl [20Tlt Kaltcne gger et al . 2011). If we 
include the possibility that giant exoplanets could har- 
bour habitable exomoons, the number of the already 
discovered potentially habitable planetary environments 
beyond the Solar System could b e rised to several tens 
(jUnderwood et al.ll2003t iKaltene ggcr 2010). Moreover, 
the existence of a plethora of other terrestrial planets 
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(TPs) and exomoons in the Galaxy is rapidly gaining 
evidence dB orucki ct al. 2011; Catan zarite fc Shaoll2011[ 
iBonfils et al.ii2011i :iKipping et al.ii2012t ). and the chances 
that a large number of potentially habitable extrasolar 
bodies could be discovered in the near future are very 
large. 

The question of which properties a planetary environ- 
ment needs in order to allow the appearance, evolution 
and diversification of life has been extens ively studied 
(for r ecent reviews see iLammer et al.l 120091 and 'Kastin^ 
|2010( ) . Two basic and complementary physical conditions 
must be fulfilled: the presence of an atmo sphere and the 
existence of liquid water on the surface (jKasting et al.l 
II993D . However, the fulfilment of these basic condi- 
tions depends on many complex and diverse endogenous 
and exogenous factors (for a comprehensive enumer a- 
tion of these factors see e.g. iWard fc Brownie^ 120001 or 
ILammer eraI1l2010t ) 

The existence and long-term stability of an intense 
planetary magnetic fi eld (PMF) is one of these rele- 
vant factors (see e.g. iGrieCmeier et all 120101 and refer- 
ences therein). It has been shown that a strong enough 
PMF could protect the atmosphere of potentially hab- 
itable planets, especially its valuable content of water 
and other volatiles, against the erosive action of the stel- 
lar wind (jLammer et al.ll2003l l2007t iKhodachenko et al.l 
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l2007t iChaufrav et al.l I2007D . Planetary magnetospheres 
also act as shields against the potentially harmful effects 
that the stellar high energy particles and galactic cos- 
mic rays (CR) produce in the life-forms evolving on the 
planetary surface (see e.g. IGrieBmeier et al.]l2005f ). Even 
in the case that life could arise and evolve on unmag- 
netized planets, the detection of atmospheric biosigna- 
tures would be also affected by a higher flux of high en- 
ergy particles including CR, especi ally if the planet is 
close to ver y active M-dwarfs (dM) (jGrenfell et al.ll2007l : 
iSegura et a l. 2010). 

It has been recently predicted that most of the TPs 
in our Galaxy coul d be fo und around dM stars (iBoss 



2006 1 iMavor fc Udrv 2008; Scab et al.ll2007HRauer eta! 
20 lit iBonfilset al.l i2011i) . Actually ~ 20% of the 



presently confirmed super-Earths belong to planetary 
systems around stars of this type. Planets in the HZ 
of low mass stars ( M^ 5, 0-QMp ) ) wou ld be tidally locked 
poshiet all 119971: iHeller et al1l20lH) . a condition that 
poses se rious limitations to their potential habitability 
(see e.g. iKite et al]|2011l and references therein). Tidally 
locked planets inside the HZ of dMs have periods in 
the range of 5 — 100 days, a condition that has com- 
monly been associated wit h the almost complete l ack of a 
protective magnetic field (jGrieBmeier et ani2004[ ). How- 
ever, the relation between rotation and PMF properties, 
that is critical at assessing the magnetic protection of 
slowly rotating planets, is more complex than previously 
thought. In particular, a detailed knowledge of the ther- 
mal evolution of the planet is required to predict not only 
the intensity but also the regime (dipolar or multipolar) 
of the PMF for a given pla netary mass and rotation rate 
(jZuluaga fc Cuartas|[2012l ). 

Although several authors have extensively studied the 
protection that intrinsic PMF would provide to extraso- 
lar p l anets (|Gr icfimcicr ct al. 2005; Khodachcnko ct all 
l2007HLam^^et al...2007; .Griefimcier et al...2009,.2010. ). 
all have disregarded the influence that thermal evolu- 
tion has in the evolution of planetary magnetic proper- 
ties. They have used also outdated dy namo scaling-laws 
that have been recently revised (see iChristensenI 120 lOl 
and references therein). The role of rotation in deter- 
mining the PMF properties that is critical in assessing 
the case of tidally locked p lanets has also been overlooked 
(jZuluaga fc Cuartas|[2012l ). 

We develop here a comprenhensive model for the evo- 
lution of the magnetic protection of potentially hab- 
itable TPs around GKM main sequence stars. To 
achieve this goal we integrate in a single framework 
a parametrized thermal evolution mod el based on 
the most recent advances i n the field (iGa idos et ^ 
120101 Tachina mi et al.l 120111 : iStamenkovic et al. 20lC 
up-to- date dy namo sc aling-laws ( ChristensenI 120101 : 
iZuluaga fc Cuarta£ll2012i) and phenomenological models 
for the evolution of the stellar wind and p lanet-star mag- 
netic interaction (jGrieCmeier et al1l2010[) . Our model is 
aimed at: 1) understanding the influence of thermal evo- 
lution in the magnetic protection of TPs, 2) assesing the 
role of low rotation periods in the evolution of the mag- 
netic protection of tidally-locked habitable planets, 3) 
placing more realistic constraints on the magnetic prop- 
erties of potentially habitable TPs suitable for future 
studies of atmospheric mass-loss or the CR effect on the 
atmospheric chemistry or on life itself and 4) estimat- 



ing by the first time the magnetic properties of already 
discovered super-Earths in the HZ of their host stars. 

Our work is a step-forward in the understanding of 
planetary magnetic protection because it puts together 
in a single model the evolution of the magnetic proper- 
ties of the planet and hence its dependence on planetary 
mass and composition and the role of the planet-star in- 
teraction into determining the resulting level of magnetic 
protection. Previous models of the former (thermal evo- 
lution and intrinsic magnetic properties) did not consider 
the interaction of the planetary magnetic field with the 
evolving stellar wind which is finally the factor that de- 
termine the actual level of planetary magnetic protection. 
On the other hand previous attempts to study the planet- 
star interaction overlooked the non-trivial dependence of 
intrinsic magnetic properties on planetary thermal evo- 
lution and hence on planetary mass and rotation rate. 
Additionally, and for the first time, we are attempting 
here to calculate the magnetic properties of the already 
discovered potentially habitable super-Earths, GJ 667Cc, 
Gl 581d and HD 40307g. 

This paper is organized as follows. Section [5] is 
aimed at introducing the properties of planetary mag- 
netospheres we should estimate in order to evaluate the 
level of magnetic protection of a potentially habitable 
terrestrial planet. Once those properties are expressed 
in terms of two basic physical quantities, the planetary 
magnetic dipole moment and the pressure of the stellar 
wind, we proceed to describe how those quantities can 
be estimated by modelling the thermal evolution of the 
planet (section ^ , scaling the dynamo properties from 
the planetary thermal and rotational properties (section 
11]) , and modelling the interaction between the star and 
the planet (section[5|). In section|6]we apply our model to 
evaluate the level of magnetic protection of hypothetical 
potentially habitable TPs as well as the already discov- 
ered habitable super-Earths. We also present there the 
results of a numerical analysis aimed at evaluating the 
sensitivity of our model to uncertainties in the compo- 
sition of the planet and to other critical parameters of 
the model. In section [7] we discuss the limitations of our 
model, present an example of the way in which our re- 
sults could be applied to estimate the mass-loss rate from 
already and yet to be discovered potentially habitable 
TPs and discuss the observational prospects to validate 
or improve the model. Finally a summary and several 
conclusions drawn from this research are presented in 
section El 

2. CRITICAL PROPERTIES OF AN EVOLVING 
MAGNETOSPHERE 

The interaction between the PMF, the interplanetary 
magnetic field (IMF) and the stellar wind creates a mag- 
netic cavity around the planet known as the magneto- 
sphere. Although magnetospheres are very complex sys- 
tems, their global properties are continuous functions of 
only two physical variables (Siscoe fc Christopher 1975): 
the magnetic dipole moment of the planet, Ai, and the 
dynamic pressure of the stellar wind, Psw Dipole mo- 
ment is defined in the multipolar expansion of the mag- 
netic field strength: 



Bpir) 
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where Bp{r) is the angular-averaged PMF strength 
measured at a distance r from the planet center and 
/xq = 477 X 10~^ H/m is the vacuum permeability. In 
the rest of the paper we will drop-off the higher order 
terms in 1/r (multipolar terms) and focus on the dipolar 
component of the field i?p'P which is explicitly given by 
the first term of the right side in eq. ([T]). 

The dynamic pressure of the stellar wind is given by 



Psw — mnv. 



off 



2nkBT. 



(2) 



where m and n are the typical mass of a wind parti- 
cle (mostly protons) and its number density, respectively. 
Here Uoff — [v^^i + ^pY^'^ i^ the effective velocity of the 
stellar wind as measured in the reference frame of the 
planet whose orbital velocity is Vp. T is the local tem- 
perature of the wind plasma and ks — 1.38 x lO^^'^ j/K 
is the Boltzmann constant. 

There arc three basic properties of planetary magne- 
tospheres we are interested in: 1) the maximum mag- 
netopause field intensity B^p, which is a proxy of the 
flux of high energy particles entering into the magneto- 
spheric cavity, 2) the standoff or stagnation radius, Rs, 
a measure of the size of the dayside magnetosphere, and 
3) the area of the polar cap Ape that measures the to- 
tal area of the planetary atmosphere exposed to open 
field lines through which particles can escape to inter- 
planetary space. The value of these quantities provides 
information about the level of exposure that a habitable 
planet has to the erosive effects of stellar wind and the 
potentially harmful effects of the CR. 

The maximum value of the magnetopause field inten- 
sity -Bmp is estimated from the balance between the mag- 
netic pressure Pmp = Bf /{2^i^) and the dynamic stellar 
wind pressure Psw (eq. ([2])), 



B„ 



(2A^o)^/^Pi/^ 



(3) 



Here we are assuming that the pressure exerted by the 
plasma inside the magnetospheric cavity is negligible (see 
discussion below). 

Although magnetopause fields arise from very complex 
processes (Chapman-Ferraro and other complex currents 
at the magnetosphere boundary), in simplified models 
-Bmp is assumed proportional to the PMF i ntensity i?D 
as measure d at the substellar point r — Rs ()Meadlll9"64 : 



Smp = 2/oPp(r = Rs) 



/oMo 

27r 



V2MRg^ 



(4) 



/o is a numerical enhancement factor of order 1 that 
can be estimated numerically. We are assuming here 
that the dipolar component of the intrinsic field (first 
term in the r.h.s. of eq. ([T])) dominates at magnetopause 
distances even in slightly dipolar PMF. 

Combining equations Q and (|4|) we estimate the 
standoff distance: 



/ f2\ 1/6 

I Aio/o \ ^i/3p_i/6 



that can be expressed in terms of the present dipole 



moment of the Earth A^© = 7.768 x 10^^ A m^ and the 
average dynamic pressure of the solar wind as m easured 
at the orbit of our planet Psw© — 2.24 x 10~^ (jStacevI 
[1991 iGriefimeier et al. 2005): 
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9.75 
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(5) 



It is important to stress that the value of Rs estimated 
with eq. ([5]) assumes a negligible value of the plasma 
pressure inside the magentospheric cavity. This approx- 
imation is valid if at least one of these conditions is ful- 
filled: 1) the planetary magnetic field is very intense, 2) 
the dynamic pressure of the stellar wind is small, or 3) 
the planetary atmosphere is not too bloated by the XUV 
radiation. In the case when any of these conditions are 
fulfilled we will refer to Rs as given by eq. ([5]) as the 
magnetic standoff distance which is an underestimation 
of the actual size of the magnetosphere. 

Last but not least we are interested in evaluating the 
area of the polar cap. This is the region in the magne- 
tosphere where magnetic field lines could be open into 
the interplanetary sp ace or to the magnetotail region. 
iSiscoe fc ChenI (|1975| ) have shown that the area of the 
polar cap Ape scales with dipole moment and dynamic 
pressure as: 



ipc 



47rP2 



= 4.63% 



M 



-1/3 
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P 



SW0 



1/6 



(6) 



Here we have normalized the polar cap area with the 
total area of the atmosphere AnRp, and assumed the at- 
mosphere has a scale-height much smaller than planetary 
radius Rp. 

In order to model the evolution of these three key mag- 
netosphere properties we need to estimate the surface 
dipolar component of the PMF Pp'P(Pp) (from which we 
can obtain the dipole moment M), the average number 
density n, velocity VeS and temperature T of the stellar 
wind (which are required to predict the dynamic pres- 
sure Psw)- These quantities depend in general on time 
and also on different planetary and stellar properties. 
In the following sections we describe our model for the 
calculation of the evolving values of these fundamental 
quantities. 

3. THERMAL AND DYNAMO EVOLUTION 

We assume here that the main source of a global PMF 
in TPs is the action of a d ynamo powered by con vection 
in a liquid metallic core ( Stevensoiil 119831 120031) . This 
assumption is reasonable since the Moon and all rocky 
planets in the solar system, regardless of their different 
origins and compositions, seem to have in the present or 
to have had in the past an iron core dynamo (see e.g. 
lStevensonll2010[ ). Other potential sources of PMFs such 
as body currents induced by the stellar magnetic field or 
dynamo action in a mantle of ice, water or magma are 
not considered here but left for future research. 

The properties and evolution of a core dynamo will de- 
pend on the internal structure and thermal history of the 
planet. Thermal evolution of TPs, specially the Earth it- 
self, has been studied for decades (for a recent review see 
iNimmol |2009() . A diversity of thermal evolution mod- 
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els for planets large r than the E arth have recently ap- 
peare d in literature (iPapuc fc Da vies 200 81 iGaidoset ahl 
2OT0t iTachinami et alJ 120111: iDriscoU fc Olson! 120111: 
Stamenkovic et al.ll2011| ) . But the lack of observational 
evidence against which we can compare the predictions 
of these models has left too much room for uncertainties 
especially regarding mantle rheology, core composition 
and thermodynamic properties. Albeit these fundamen- 
tal limitations, a global picture of the thermal history 
of super-Earths has started to arise. Here we follow the 
lines of Labrossc 2003 and Gaidos et al. 2010 develop- 
ing a parametrized thermal evolution model which com- 
bines a simplified model of the interior structure and a 
parametrized description of the core and mantle rheol- 
ogy- 

Our model includes several distinctive characteristics 
in comparison to previous ones. The most important one 
is an up-to-date treatment of mantle rheology. For that 
purpose we use two different formulae to compute the vis- 
cosity of the upper and lower mantle. By lower mantle 
we understand here the region of the mantle close to the 
core mantle boundary (CMB). This is the hottest part 
of the mantle. The upper mantle is the outer cold part 
of this layer. It is customary to describe both regions 
with the same rheology albeit their very different miner- 
alogical compositions. Additionally thermal and density 
profiles in the mantle following the same prescription as 
in the core. We also use a different ansatz to assign initial 
values to lower mantle temperature and to the tempera- 
ture contrast across the CMB, two of the most uncertain 
quantities in thermal evolution models. Using our ansatz 
we avoid assigning arbitrary initial values to these criti- 
cal parameters but more importantly we are able to find 
a unified method to set the value of these temperatures 
in planets with very different masses. It is also impor- 
tant to notice that in other models these temperatures 
were set by hand or were treated as free parameters in 
the model 

Four key properties should be predicted by any ther- 
mal evolution model in order to calculate the magnetic 
properties of a planet: 1) the total available convective 
power Qconv, providing the energy required for magnetic 
field amplification through dynamo action; 2) the radius 
of the solid inner core i?ic and from there the height 
D Ri Re — Ric of the convecting shell where the dynamo 
action takes place {Re is the radius of the core) ; 3) the 
time of inner core formation tic and 4) the total dynamo 
life-time tdyn- 

In order to calculate these quantities we solve sim- 
ply parametrized energy and entropy equations of 
balance describing the flux of heat and entropy in 
the planetary core and mantle. As stated before 
our model is based on the interior structure model 
by VAL06) and in thermal evolution ni o dels previ- 
ously developed by [Schubert et al.l IT9791 iSteyensonI 



1983, 'Nimmo fc Stevenson' '2000*. 'Labrosse et al 



Labrossc 2003, Gubbins et al. 2003, Gub bins et al.l 



20011 
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20 111 For a detailed description of the fundamental 
physics behind the thermal evolution model developed 
here please refer to these earlier studies. 

3.1. Interior structure 



Our one-dimensional model for the interior assumes 
a planet made by two well differentiated chemically 
and mineralogically homogeneous shells: a rocky man- 
tle made out of olivine and perovskite and a core made 
by iron plus other light elements. 

The mechanical conditions inside the planet (pres- 
sure P, density p and gravitational field g) are com- 
puted by solving simultaneously the continuity, Adams- 
Williamson and hydrostatic equillibrium equations (eqs. 
(l)-(4) in VAL06). For all planetary masses we assume 
boundary conditions, p{r = Rp) = 4000 kg m~'^ and 
P{r = Rp) = Pa. For each planetary mass and core 
mass fraction CMF= Mcorc/M^p we use a RK4 integrator 
and a shooting method to compute consistently the core 
Re and planetary radius Rp. 

For the sake of simplicity, we do not include in the in- 
terior model the two- or even three-layered structure of 
the mantle. Instead of that we assume a mantle com- 
pletely made of perovskite-postperovskite (ppv). This is 
the reason why we take p(r = Rp) = 4000 kg m~'^ in- 
stead of the more realistic value of 3000 kg m""^. With 
a single layer and a realistic surface density our model is 
able to reproduce the present interior properties of the 
Earth. 

In all cases we use the Vinet equation of state instead of 
the commonly used third order Birch-Murnaghan equa- 
tion (BM3). It is well known that the BM3 follows from 
a finite strain expansion and does not accurately pre- 
dict the properties of the material for the typical pres- 
sures found in super-Earths , i. e. 100 — 1000 GPa fo r 
Mp = 1 - lOMe (VAL06 and ITachinami et al.l [20Tll) . 
We have ignored thermal corrections to the adiabatic 
compressibility, i.e. Ks{p, T) w Ks{p, 300 K) + AK,{T) 
(VAL06). This assumption allows us to decouple at run- 
time the CPU intensive calculation of the thermal profile 
from the mechanical structure at each time-step in the 
thermal evolution integration. 

Although we have ignored the "first-order" effect of the 
temperature in the mechanical structure, we have taken 
into account "second-order" effects produced by phase 
transitions inside the mantle and core. Using the ini- 
tial temperature profile inside the mantle we calculate 
the radius of transition from olivine to perovskite (ne- 
glecting the effect of an intermediate layer of wadsleyite). 
For that purpose we use the reduced pressure function 11 



jrpc 

i3i 



(|Christensenlll985l : rWernsteinlll992tlValencia et al.ll2007| ). 
For simplicity the position of the transition layer is as- 
sumed constant during the whole thermal evolution of 
the planet. We have verified that this assumption does 
not change significantly the mechanical properties inside 
the planet at least not at a level affecting the thermal 
evolution itself. 

Inside the core we update continuously the radius of 
the transition from solid to liquid iron (see below). For 
that purpose we use the thermal profile computed at the 
previous time-step. To avoid a continuous update of the 
mechanical structure we assume in all the cases that at 
the transition from solid to liquid the density of iron 
changes by a constant factor Ap = (ps — pi)/pi- Here 
the reference density of the solid ps (when applied) is 
computed using the Vinet equation evaluated at a point 
in the very center of the planet. 

Table [1] enumerates the relevant physical parameters 
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of the interior structure and thermal evolution model. 

3.2. Core thermal evolution 

In order to compute the thermal evolution of the core 
we solve the e quatio ns of en ergy and entropy balance 
(jLabrosse et al.ll200UlNimmoll2009D : 



where Dc = ^j2>Cp/2TTapcG is the scale height of tem- 
perature, a is the isothermal expansivity (assumed for 
simplicity constant along the core) and pc is the density 
at core center. Using this fit the bulk secular heat ca- 
pacity Cs = Qs/{Mc dTc/dt) can be obtained from eq. 



Qc = Q. + MQg + Qi) 



^ = Ei+f,{Es+Eg)~Ek. 



(7) 



(8) 



Here Qc is the total heat flowing through the CMB and 
$ is the total entropy dissipated in the core. Eg and Qs 
are the entropy and heat released by the secular cooling, 
Eg and Qg are the contribution to entropy and heat due 
to the redistribution of gravitational potential when light 
elements are released at the liquid-solid interface (buoy- 
ant energy) , Ei and Qi are the entropy and heat released 
by the phase transition (latent heat) and E]^ is a term 
accounting for the sink of entropy due to the conduction 
of heat along the core. We have avoided the terms com- 
ing from radioactive and pressure heating because their 
contribution are negligibl e at t he typical conditions in- 
side super-Earths (Nimmoll2009[ ). As long as the buoyant 
and latent entropy and heat are only present when a solid 
inner core exists we have introduced a boolean variable 
fi that turn-on these terms when the condition for the 
solidification of the inner core arises. 

The terms in the energy and entropy balance are a 
function of the time-derivative of the temperature profile 
dT(r,t)/ dt (for detaile d expressions of these terms see 
table 1 in lNimmoll2009| ) . As an example the secular heat 
and entropy are given by: 



Qs 

E_. 



pcp — dV, 



pcp 



dt 
1 

Tjt) T{r,t) 



1 



dT{r,t) 
dt 



dV. (9) 



Here Cp is the specific heat of the core alloy and T^ 
is the temperature at the CMB. If we assume that the 
temperature profile of the core does not change dur- 
ing the thermal evolution, we can write temperature as 
T{r,t) — fc{r)Tc{t). Here fc{r) is the core temperature 
radial profile that we will assume adiabatic (see below). 
It should be noted again that Tc{t) — T{r — Rc,t). 

With this assumption the energy balance equation ^ 
can be written as a first order differential equation on Tc'. 

Jrr\ 

Q, - M, [Cs + h {Cg + Ci)]^ (10) 

where Mc is the total mass of the core and Cg, Cg and 
Ci are core bulk heat capacities which can be expressed 
as volumetric integrals of the radial profile fc{r). In this 
equation the total heat Qc is intrinsically a function of 
Tc and should be computed independently (see below). 

Using a si mple exponentia l fit fo r the core density, as 
proposed by iLabrosse et al.l ()200lD , the adiabatic tern - 
perature profile can be approximated as (|Labrossell2003D : 



fair) = exp 
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(11) 
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r dr 



(12) 



Analogous expressions for Cg and Ci are obtai ned from 
the de finition of Qg and Qi as given in table 1 of lNimmol 
((200a . 

The total heat released by the core Qc{Tc) is calcu- 
lated herejasingjhe boundary layer theory (BLT) (see 
e.g. IStevensonlll983[). Under this approximation Qc is 
given by (|Ricardll2009( ) : 



Qc = 47ri?cA:,„ArcMBNuc, 



(13) 



where km is the thermal conductivity of the lower 
mantle, ATcmb = Tc — Ti is the temperature contrast 
across the CMB, T; is the lower mantle temperature, and 
Nuc ^ (Rae/Ra^,)^/^ is the Nusselt number at the core 
(jSchubert et al.|[2001l ) . The critical Rayleigh number Ra* 
is a free parameter in our model (see table [ij . 

The local Rayleigh number Rac at the CMB is calcu- 
lated under the ass umption of a boundary heated from 
below (|Ricardll2009t) . 



Rac 



p g a ATcMBJRp - Rcf 



(14) 



where g is the gravitational field and Kc the thermal 
diffusivity at the CMB. The value of the dynamic viscos- 
ity 77c, which is strongly dependent on temperature, could 
be suita bly computed using the so-called film tempera- 
ture (see lManga et al.ll2001l and reference therein). This 
temperature could be computed in general as a weighted 
average of the temperatures at the boundaries. 



T,,c^icTc + {l-£.c)Tu 



(15) 



where the weighting coefficient S,c is a free parameter 
whose value is chosen in order to reproduce the thermal 
properties of the Earth (see table [IJ. 

To model the formation and evolution of the solid inner 
core we need to compare at each time the temperature 
profile with the iron solidus. We use here the Lindemann 
law as parametrized by VAL06: 



dlogT 
dlogp 



2 [7 - S{p)] 



(16) 



Here 6{p) w 1/3 and 7 is an effective Griineisen pa- 
rameter that is assumed for simplicity constant. To inte- 
grate this equation we use the numerical density profile 
provided by the interior model and the reference values 
po — 8300 kg/m (pure iron) and tq = 1808 K. 

The central temperature T{r — 0, t) and the solidus 
at that point T{r = 0) are compared at each time step. 
When T(0,tic) w r(0) {tic is the time of inner core for- 
mation) we turn-on the buoyant and latent heat terms 
in eqs. ^ and ^, i.e. set fi = 1, and continue the in- 
tegration including these terms. The radius of the inner 
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core at times t > tjc is obtaine d by solving the equa- 
tio n proposed by iNimmol ()2009l ) and further developed 
bv lGaidos erair(|2010n . 



IT 



Dl 



1 dTr 



2i?ic(A- l)Tc dt 



(17) 



Here A is the ratio between the gradient of the solidus 
(eq. ([16])) and the actual temperature gradient Tc{t)fc{r) 
as measured at -Ric(i)- 

When the core cools down below a given level the outer 
layers start to stratificate. Here we model the effect of 
stratification by correcting the radius and temperature 
of the core following the prescription by iGaidos et alj 
()2010[) . When stratified the eff ective radius of the core 
is reduced to Ri, (eq. (27) in IGaidos et al.l I2010D and 
the te mperature at the co re surface is increase to T^, (eq. 
(28) in lGaidos et al.|[2010D . The stratification of the core 
reduces the height of the convective shell which leads to a 
reduction of the Coriolis force potentially enhancing the 
intensity of the dynamo-generated magnetic field. 

The estimation of the dynamo properties requires the 
computation of the available convective power Qconv 
Qconv is calculated here assuming that most of the dissi- 
pation occurs at the top of the core. Under this assump- 
tion, 



,{t) « <i>{t)T,{t) 



(18) 



where the total entropy $ is computed from the En- 
tropy balance (eq. ([5])) using the solution for the tem- 
perature profile Tc{t) fc{r) . 

When $(i) becomes negative, i.e. Ei + Eg + Eg < Ek 
in eq. ([8|), Qconv gets also negative and convection is 
not longer efficient to transport energy across the outer 
core. Under this condition the dynamo is shut down. 
The integration stops when this condition is fulfilled at 
a time we label as the dynamo life-time tdyn- 

3.3. Mantle thermal evolution 

One of the novel features of our thermal evolution 
model is that we treat mantle thermal evolution with a 
similar formalism as that described before for the metal- 
lic core. 

The energy balance in the mantle can be written as: 



Qrn = XrQr +Qs+Qc 



(19) 



Here Qm is the total heat fiowing out through the sur- 
face boundary (SB), Qr is the heat produced in the decay 
of radioactive nuclides inside the mantle, Qs is the sec- 
ular heat and Qc is the heat coming from the core (eq. 

(USD). 

We use here the standard expressions and parame- 
ters for the radi oactive energy production as given by 
iKite et all ()2009() . However, in order to correct for the 
non- homogeneous distribution of radioactive elements in 
the mantle, we introduce a multiplicative correction fac- 
tor Xr- Here we adopt Xr = 1.253 that fits well the Earth 
properties. We have verified that the thermal evolution 
is not too sensitive to Xr and have assumed the same 
value for all planetary masses. 

The secular heat in the mantle is computed using an 
analogous expression to eq. (HJ. As in the case of the 
core we assume that the temperature radial profile does 



not change during the thermal evolution. Under this 
assumption the temperature profile in the mantle can 
be also written as T{r,t) = Tm{t) fmir) ■ In this case 
Tm{t) = T{r = Rp,t) is the temperature right below the 
surface boundary layer (see figure [T]). 

Assuming an adiabatical temperature profile in the 
mantle we can also write: 



fm{r) = exp 



Rl 



£)2 



in analogy to the core temperature profile, eq. (jlip . 
In this expression Dm is the temperature scale height for 
the mantle which is re lated to the dens ity scale height Lm 
through D^ = L^ /^ [Labrossd (|2003f ). In our simplified 
model we take the values of the density at the boundaries 
of the mantle and obtain analytically an estimate for Lm 
and hence for Dm- 

The energy balance in the mantle is balanced when 
we independently calculate the heat Q,„ at the SB as 
a function of Tm- In this case the presence or not of 
mobile lids play an important role into determining the 
efficiency with which the planet gets rid of the heat com- 
ing from the mantle. In the mobile lid regime we assume 
that the outer layer is fully convective and use the BLT 
approximation to calculate Qmi 



Q. 



ML 



47ri?gfc^AT„Nu, 



(20) 



where AT^ — T„i — Ts is the temperature contrast 
across the SB and Tg the surface temperature. Since we 
are studying the thermal evolution of habitable planets 
we assume in all cases T, = 290 K. Planetary interior 
structure and thermal evolution are not too sensitive to 
surface temperature. We have verified that results are 
nearly the same for surface temperature in the range of 
250—370 K. In the mobile lid regime Nu™ obeys the same 
relationship with the critical Rayleigh number as in the 
core. In this case however we compute the local Rayleigh 
number under the assumption of material heated from 



inside (jGaidos et al.l 



assunrp 

][2oii. 



Ra„ 



igp^HjRp - R,f 



(21) 



where H = {Qr + Qc)/Mm is the density of heat in- 
side the mantle, k^ and k^ are the thermal conductivity 
and diffusivity respectively and rjm is the upper mantle 
viscosity. 

In the stagnant lid regime the SB provides a rigid 
boundary for the heat flux. In this case we adopt the 
approximation used by Nimmo & Stevens on., 2000, : 



Q 



SL 



^ttRI^ 






1/3 



p-4/3 



(22) 



Here T = —dlnrjm/dTm measures the viscosity depen- 
dence on temperature evaluated at the average mantle 
pressure. 

With all this elements at hand the energy balance at 
eq. (J19p is finally transformed into an ordinary differen- 
tial equation for the upper mantle temperature Tm(t), 



The influence of thermal evolution in the magnetic protection of terrestrial planets 



"^ \r^r 1 ^c I ^r: 



' dt 



(23) 



where Cm is the bulk heat capacity of the upper mantle 
which is calculated with an analogous expression to eq. 

m 

3.4. Initial conditions 
In order to solve the coupled differential eqs. (fTO|) . 



(|T7j) and p3|) we need to choose a proper set of initial 
conditions. 

The initial value of the upper mantle temperature 
is cho sen using the prescription by iStamenkovic et al.l 
(|2011| ). According to this prescription Tm{t = 0) is com- 
puted by integrating the pressure-dependent adiabatic 
equation up to the average pressure inside the mantle 
(Pm), 



Tm{t = 0) = 0exp 



•{^™> 



70 



Ks{P') 



dP' 



(24) 



Here 9 = 1700 K is a potential temperature 
which is assumed the s ame for all planetary masses 
(jStamenkovic et aLll201lD . Using Tm (t = 0) and the adi- 
abatic temperature profile we can get the initial lower 
mantle temperature Ti(t = 0). 

The initial value of core temperature Tc(t — 0) is 
one of the most uncertain parameters in thermal evo- 
lution models. Although nobody knows its actual value 
or its dependence on the formation history and plane- 
tary mass, it is reasonable to start the integration of a 
simplified thermal evolution model when the core tem- 
perature is of the same order as the melting point for 
MgSiOa at the lower mantle pressure. A small arbi- 
trary temperature contrast against this r eference value 
(jGaidos et al.ll2010l : lTachinami et al.ll201lD or more com- 



plicat ed mass-dependent assumptions (jPapuc fc DaviesI 
1200 8!) have been used in previous models to set the ini- 
tial core temperature. We use here a simple prescription 
that agrees reasonably well with previous attempts and 
provides a unified expression that could be used consis- 
tently for all planetary masses. 

According to our prescription the temperature contrast 
across the CMB is assumed proportional to the temper- 
ature contrast across the whole mantle, i.e. ATcmb — 
CadbATadb = £a.dh(Tm — Ti). Wc havc found that the ther- 
mal evolution properties of Earth are reproduced when 
we set eadb = 0.7. 

Using this prescription the initial core temperature is 
finally calculated using: 



T,{t - 0) = T( + fadbAT, 



adb 



(25) 



We have observed that the value of the Tc{t — 0) ob- 
tained with this prescription is very close to the per- 
ovskite melting temperature at the CMB for all the plan- 
etary masses studied here. This result shows that al- 
though our criterium is not particularly better physically 
rooted than those used in previous models, it still relies 
in just one free parameter, i.e. the ratio of mantle and 
CMB contrasts eadb- 



3.5. Rheological model 

One of the most controversial aspects and probably 
the largest source of uncertainties in thermal evolution 
models is the calculation of the rheological properties of 
sillicates and iron at high pressures and temperatures. A 
detailed discussion on this important topic is out of the 
scope of this paper. An up-to-date discussion and analy- 
sis of the dependence on pressure and temperature of vis- 
cosity in super-Earths and its infiuence in thermal evolu- 
tion c an be found in the recent wor ks bv lTachinami et al.l 
120111 and IStamenkovic et al.|[20Tll 

We use here two different models to calculate vis- 
cosity under different ranges of temperatures and pres- 
sures. For the high pressures and temperatures of 
the lower mantle we us e a Nabarro-Herring model 
(jYamazaki fc Karatol[200l . 



VMP^T) 



Rgd" 



DqAitLj, 



-rp(P,T)exp 



Here Ra 



8.31 JmoU^K"^ 



T 

(26) 
is the gas constant, d is 
the grain size and m the growing exponent, A and b 
are free parameters, Dq is the pre-exponential diffusion 
coefficient and mmoi is the molar density of perovskite 
and Tmoit(-P) is the melting temperature of perovskite 
that can be computed with the empirical fit: 



T^.\t{P) 



1=0 



All the parameters used in the visocisty model, in- 
cluding the expansion coefficients a^ in the melting tem- 
perature formula, were t aken from the recent work by 
IStamenkovic et al] (|2011f ) . The Nabarro-Herring formula 
allows us to compute rjc = T/jsmlTric) , where T^/c is the film 
temperature computed using the average in eq. (|15p . 

The upper mantle has a completely different mineral- 
ogy and it is under the influence of lower pressures and 
temperatures. Although previous works have used the 
same model and p arameters to calculate vis cosity across 
the whole mantle ([Stamenkovic et ahl 120111 for example 
use the perovskite viscosity parameters also in the olivine 
upper mantle), we have found here that using a different 
rheological model in the upper and lower mantle avoids 
an under and overestimation, respectively, of the value 
of viscosity that could have a significant effect on the 
thermal evolution. 

In the upper mantle we find that using an Arrhenius- 
type model leads to better estimates of viscosity than 
that obtained using the Nabarro-Herring model. In the 
upper mantle the Nabarro-Herring formula (which is best 
suited to describe the dependence on viscosity at high- 
pressures and temperatures) leads to huge underestima- 
tions of viscosity in that region. In the case of the Earth 
this underestimation produces values of the total mantle 
too high as compared to that observed in our planet mak- 
ing impossible to fit the thermal evolution of a simulated 
Earth. 

For the Arrhenius-type fo rmula we use th e same 
parametrization given by .Tachinami et all (|201in : 
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Va{P,T) 



51/n 



exp 



E* + PV* 
nRgT 



;(!-")/« 



(27) 



ROi — CroI Pc 



- -l/6p-2/3r)-l/3T/-l/2 



'R 



D' 



'V- 



'qV^^^P'/'. (29) 



where e is the strain rate, n is the creep index, B is 
the Barger coefficient, and E* and V* are the activa- 
tion energy and volume. The values assumed here for 
these parameters are_the same as that given in table 4 of 
[^chinami ct al. 2011| except for the activation volume 
whose value we assume here V* = 2.5 x 10~^ m^ mol~^. 
Using the formula in eq. (I27[) . the upper mantle viscosity 
is computed as r/^ = ?7A((-Rm),T'm)- 

A summary of the parameters used by our interior and 
thermal evolution models are presented in table [TJ The 
values listed in column 3 define what we will call the ref- 
erence thermal evolution model (RTEM). These reference 
values have been mostly obtained by fitting the present 
interior properties of the Earth and the global features of 
its thermal, dynamo and magnetic field evolution (time 
of inner core formation and present values of Ric , Qm 
and surface magnetic field intensit y). For the s tagnant 
lid case we use as suggested by Gaid os et al.ll2010l the val- 
ues of the parameters that globally reproduce the present 
thermal and magnetic properties of Venus. 

Figure [2] shows the results of applying the RTEM to a 
set of hypothetical TPs in the mass-range Mp = 0.5 — 
6Me. 

4. PLANETARY MAGNETIC FIELD 

In recent years improved numerical experiments have 
constrained the full set of possible scaling laws used to 
predict the properties of planetary and ste llar convection- 
driven dynamos (see iChristensenI 120101 and references 
therein). It has been found that in a wide range of 
physical conditions the global properties of a planetary 
dynamo can be expressed in terms of simple power-law 
functions of the total convective power Qconv and the size 
of the convective region. 

One of the most important results of power-based 
scaling laws is the fact that the volume averaged mag- 
netic field intensity B^^-^^ = (l/V) J B^dV does not de- 
pend on the rotation ra te of the planet (cq. (6) in 
IZuluaga fc Cuartas|[20T2l ). 



Br, 



C^rr..t^,'%''\D/Vf/'Qlil 



(28) 



Here Cerms is a fitting constant obtained from numer- 
ical dynamo experiments and its value is different in the 
case of dipolar dominated dynamos, Cg'^.^^^^ == 0.24, and 
multipolar dynamos, C^^^s = 0.18. p^, D = R^ — Ric 
and V — 47r(_R^ — Rfc)/3 are the average density, height 
and volume of the convective shell. 

The dipolar field intensity at the planetary surface, 
and hence the dipole moment of the PMF, can be esti- 
mated if we have information about the power spectrum 
of the magnetic field at the core surface. Although we 
cannot predict the relative contribution of each mode to 
the total core field strength, numerical dynamos exhibit 
an interesting property: there is a scalable dimension- 
less quantity, the local Rossby number Ro"^ , that could 
be used to distinguish dipolar dominated from multipo- 
lar dynamos. The scaling relation for Ro* is (eq. (5) in 
[Zuluaga fc Cuartas, 2012, ): 



Here CroI = 0.67 is a fitting constant and P is the 
period of rotation. It has been found that dipolar dom- 
inated fields arise systematically when dynamos have 
Ro* < 0.1. Multipolar fields arise in dynamos with val- 
ues of the local Rossby number close to and larger than 
this critical value. From eq. (|29p we see that in gen- 
eral fast rotating dynamos (low P) have dipolar dom- 
inated core fields while slowly rotating ones (large P) 
produce multipolar fields and hence fields with a much 
lower dipole moment. 

It is important to stress that the almost independence 
of Brms on rotation rate, together with the role that rota- 
tion has in the determination of the core field regime, im- 
plies that even very slowly rotating planets could have a 
magnetic energy budget of comparable sized than rapidly 
rotating planets with similar size and thermal histories. 
In the former case the magnetic energy will be redis- 
tributed among other multipolar modes rendering the 
core field more complex in space and probably also in 
time. Together all these facts introduce a non-trivial 
dependence of dipole moment on rotation rate very dif- 
ferent from that obtained with the t raditional scaling 
laws use d in previous works (see e.g. iGriefimeier et al.1 
120041 and lKhodachenko et"aIll2007D . Here we want to em- 
phasize a property that was also previously overlooked. 
Multipolar dominated dynamos produces magnetic fields 
that decays more rapidly with distance than dipolar fields 
and so it is expected that a planet with a multipolar 
magnetic field will be less protected than those having 
strongly dipolar fields. 

Using the value of Brms and Roi we can compute the 
maximum dipolar component of the field at core surface. 
For this purpose we use an upper bound to the dipolar- 
ity fraction /dip (the ratio of the dipolar component to 
the total field strength at core surface). Dipolar dom- 
inated dynamos have by definition /dip < f^^p^ — 1-0. 
The case of reversing dipolar and multipolar dynamos 
is more complex. Numerical dynamo experiments show 
that multipolar dynamos have Ro^ ^0.1 and /dip < 
fdin^ ~ 0.35. However to avoid inhomogeneities in the 
transition region around i?0;* « 0.1 we calculate a max- 
imum dipolarity fraction through a "soft step function" , 
/™°^ = a + ;3/{exp[(i?0;* - 0.l)/6] + 1} with a, /3 and 
S numerical constants that fits the envelop of the nu- 
merical dynamo data (s ee upper panel of figure 1 in 
IZuluaga fc Cuartasll20ia ). 

To connect this ratio to the volumetric averaged mag- 
netic field i?rms we use the volumetric dipolarity frac- 
tion 6dip that it is found, as shown by numerical exper- 
iments, conveniently rel ated with the maxium v alue of 
/dip through eq. (12) in IZuluaga fc Cuartasll20l3 . 



rrnin rmaa 

"dip ~ CbdipJdip 



-11/10 



(30) 



where Cbdip ~ 2.5 is again a fitting constant. It 
is important to notice here that the exponent 11/10 
is the ratio of the smallest integers close to the nu- 
merical value of the fitt ing exponent (see figure 1 in 
IZuluaga fc C uartas '2012^. We use this convention fol- 
lowing lOlson fc Christcnscn.,2006, . 
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Parameter 






Definition 


Value 




Ref. 


Bulk 


CMF 

Ps 








Core mass fraction 
Surface temperature 
Surface pressure 


0.325 
290 K 
bar 




- 


Inner core 


Po, Ko, 

fee 

AS 


^0. 


70, 


9, So 


Material 

Equation of state parameters 
Thermal conductivity 
Entropy of fusion 


Fe 

8300 kg m"'', 160.2 GPa, 5.82, 

40 Wm-i K-i 

118 j kg-^K-i 


1.36, 0.91, 998 K 


A 
A 
B 
C 


Outer core 



- 






Material 


FC(o,8)I^<^S(o^2) 




A 


Po, Ko, Kg, 70, 


q, 


Bo 


Equation of state parameters 


7171 kg m-^, 150.2 GPa, 


5.675, 1.36, 0.91, 998 K 


A 


a 






Thermal expansivity 


1.4 xlO"^ K-^ 




D 


Cp 






Specific heat 


850 j kg-^ K-i 




C 


fee 






Thermal conductivity 


40 W m"^ K"^ 




B 


«c 






Thermal diffusivity 


6.5 X 10"^ m^ s-i 




E 


AS 






Entropy of fusion 


118 j kg-i K-^ 




C 


<^a.db 






Adiabatic factor for Tc{t ^ 0) 


0.7 




- 


£c 






Weight of Tc in core viscosity 


0.4 




- 


Lower mantle 


- 






Material 


pv+fniw 




A 


Po, Ko, Ko, 70, 


q. 


Bo 


Equation of state parameters 


4152 kg m--"*, 223.6 GPa, 


4.274, 1.48, 1.4, 1070 K 


A 


d, m. A, b. Do, 


m 


mol 


Viscosity parameters 


1x10"^ m, 2, 13.3, 12.33 


, 2.7x10"^° m^ s^S 0.10039 kg mol^^ 


F 


a 






Thermal expansivity 


2.4 xlO"'^K"^ 




D 


Cp 






Specific heat 


1250 jkg-^ K-i 




C 


Km 






Thermal conductivity 


6 W m-^ K"^ 




C 


f^m 






Thermal diffusivity 


7.5 X 10~^ m^s"^ 




E 


AS 






Entropy of fusion 


130 j kg-^ K-i 




C 


Upper mantle 


- 






Material 


olivine 




A 


Po, Ko, Kg, 70, 


q. 


Bo 


Equation of state parameters 


3347 kg m"^, 126.8 GPa, 


4.274, 0.99, 2.1, 809 K 


A 


B, n, E*, e 






Viscosity parameters 


3.5 X 10"^^ Pa""s"\ 3, 


430x10^ j mol-\ 1 X 10"^^ s-^ 


D 


V 






Activation volume 


2.5 xlO"'^ m^ mor^ 




F 


a 






Thermal expansivity 


3.6 xlO"'^ K"^ 




D 


Cp 






Specific heat 


1250 j kg-i K-i 




C 


Km 






Thermal conductivity 


6 W m"^ K-^ 




C 


Km 






Thermal diffusivity 


7.5 X 10"^ m^ s"i 




E 


e 






Potential temperature 


1700 K 




F 


Xr 






Radiactive heat correction 


1.253 




- 



TABLE 1 
Reference Thermal E volution Model parame ters. Sources: (A) VAL 06. fBI INimmoI (120091 ) . CCiIGaidos et al] H2010() . (D) 
ITachinami et AL.I (120111) . fEI IRicardI 1120091) . (F) IStamenkovic et alT 



~ mrv[u 



Finally by combining eqs. ([28|) - p0p we can compute 
an upper bound to the dipolar component of the field at 
the CMB: 



1 J^rna 

Bf^P < iSrms = ^ 



rnaxll/lO 



^dip 



By 



Cbdii 



(31) 



The surface dipolar field strength is estimated using, 



Bf^R,) = BfP ( 1^ 



(32) 



and finally the total dipole moment is calculated using 
eq. dl]) for r = Rp. 

It should be emphasized that the surface magnetic field 
intensity determined using eq. (J32]) overestimates the 
PMF dipolar component. The actual field could be much 
more complex spatially and the dipolar component could 
be lower. As a consequence our model can only pre- 
dict the maximum level of protection that a given planet 
could have from a dynamo-generated intrinsic PMF. 

The results of applying the RTEM to calculate the 
properties of the magnetic field of TPs in the mass range 
0.5-4.0 M® using the scaling laws in eqs. (gH), (gS]) and 
(|3ip are summarized in figures |3] and 21 In figure |3] we 
show the local Rossby number, the maximum dipolar 
field intensity and the dipole moment as a function of 



time computed for planets with different mass and two 
different periods of rotation {P — 1 day and P = 2 days). 
This figure shows the effect that rotation has on the evo- 
lution of dynamo geometry and hence in the maximum 
attainable dipolar field intensity at the planetary surface. 
In figure [5] we have summarized in mass-period (M-P) 
diagrams (Zuluaga & Cuartas 20i3) the evolution of the 
dipole moment for planets with long-lived dynamos. We 
see that for periods lower than 1 day and larger than 5-7 
days the dipole moment is nearly independent of rota- 
tion. Slowly rotating planets have a non-negligible dipole 
moment which is systematically larger for more massive 
planets. 

5. PLANET-STAR INTERACTION 

The PMF properties constrained using the thermal 
evolution model and the dynamo-scaling laws are not 
enough to evaluate the level of magnetic protection of a 
potentially habitable TP. We also need to estimate also 
the magnetosphere and stellar properties (stellar wind 
and luminosity) as a function of time in order to assess 
properly the level of star-planet interaction. 

Since the model developed in previous sections pro- 
vides only the maximum intensity of the PMF, we will 
be interested here into constraint the magnetopshere and 
stellar properties from below, i.e. to find the lower level 
of "stellar aggression" for a given star-planet configura- 
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tion. Combining upper bounds of PMF properties and 
lower bounds for the star-planet interaction will produce 
an overestimation of the overall magnetic protection of a 
planet. If under this model a given star-planet configu- 
ration is not suitable to provide enough magnetic protec- 
tion to the planet, the actual case should be much worse. 
But if, on the other hand, our upper-limit approach pre- 
dicts a high level of magnetic protection, the actual case 
could still be that of an unprotected planet. Therefore 
our model is capable at predicting which planets will be 
unprotected but less able when predicting which ones will 
be actually protected. 

5.1. The Habitable Zone (HZ) and tidally locking limits 

Surface temperature and hence "first-order" habitabil- 
ity of a planet depends on three basic factors: 1) the 
fundamental properties of the star (luminosity Li,, effec- 
tive temperature T^, and radius Ri,) 2) the average star- 
planet distance (distance to the HZ) and 3) the conmen- 
surability of planetary rotation and orbital period (tidal 
locking). These properties should be properly modelled 
in order to assess the degree of star-planet interaction 
which are critical at determining the magnetic protec- 
tion. 

The basic properties of main-sequence stars of different 
masses and metallicities have been studied for decades 
and are becoming critical in assessing the actual proper- 
ties of newly discovered exoplanets. The case of low mass 
main sequence stars (GKM) are particularly important 
in providing the properties of the stars with the highest 
potential to harbor habitable planets with evolved and 
diverse biospheres. 

In this work we will use the theoretical results by 
IBaraffe et all 119981 (hereafter BAR98) that predict the 
evolution of different metallicities main-sequence GKM 
stars. We have chosen from that model those results cor- 
responding to the case of solar metallicity stars. We have 
disregarded the fact that the basic stellar properties ac- 
tually evolve during the critical period where magnetic 
protection will be evaluated, i.e. i = 0.5 — 3 Gyr. To 
be consistent with the purpose of estimating upper lim- 
its to magnetic protection, we took the stellar proper- 
ties as provided by the model at the highest end of the 
time interval, i.e. t = 3 Gyr. Since luminosity increases 
with time in GKM stars this assumption guarantees the 
largest distance of the HZ and hence the lowest effects of 
the stellar insolation and the stellar wind. 

In order to estimate the HZ limits we use the r ecentl y 
updated values calculated by IKopparapu et al.l (|2013D . 
In particular we use the interpolation formula in eq. (2) 
and coefficients in Table (2) to compute the most con- 
servative limits of Recent Venus and Early Mars. The 
limits calculated for the stellar properties assumed here 
are depicted in figure [S] 

The orbital and rotational properties of planets at 
close-in orbits are strongly affected by the gravitational 
and tidal interaction with the host star. Tidal torque 
dampens the primordial rotation and axis tilt leaving the 
planet in a final resonant equilibrium where the period 
of rotation P becomes commensurable with the orbital 
period Pq, 



Here n is an integer larger than or equal to 2. The 
value of n is determined by multiple dynamic factors, 
the most importa nt being the orbital planetary ec- 
centricity fc econt e et al.l 120101: iFerraz-Mello et al.|[200l 
IHeller et al.l 120111) . In the solar system the tidal inter- 
action between the Sun and Mercury has trapped the 
planet in a 3:2 resonance. In the case of GL 581d, de- 
tailed dynamic models p redict a resonant 2:1 equilibrium 
state (jHeller et al.l 1201 ll ). i.e. the rotation period of the 
planet is a half of its orbital period. 

Although estimating in general the time required for 
the "tidal erosion" is very hard given the large uncer- 
tainties in the k ey physical parameters involved (see 
IHeller et al.ll201J] for a detailed discussion) the maximum 
distance Otid at which a solid planet in a circular orbit be- 
comes tidally locked before a given time t can be roughly 
estimated by (jPealelll977l) : 



atid(<) = 0.5AU 



(M,/MQ)2Pp,im 

Q 



1/6 



-1/6 



(34) 



Here the primordial period of rotation Pprim 



P:Po 



(33) 



should be 
expressed in hours, t in Gyr and Q is the dimensionless 
dissipation function. For the purposes of this work we 
assume a primordial period of rotatio n Pprim — 17 hours 
(jVarga et al.lll998t iDenis et"all 120111 ) and a dissipation 
funct ion Q « 100 (jHenning et al.l [20091 IHeller et al.l 
"Ml!). 

In figure [5] we summarize the properties of solar metal- 
licity GKM main sequence stars provided by the BAR98 
model and the corresponding limits of the HZ and tidal 
locking maximum distance. The properties of the host 
stars of the already discovered potentially habitable 
super-Earths, GL 581d, GJ 667Cc and HD 40307g, are 
also highlighted in this figure. 

5.2. Stellar wind 

The stellar wind and cosmic rays pose the highest risks 
for a magnetically unprotected potentially habitable ter- 
restrial planet. The dynamic pressure of the wind is able 
to obliterate an exposed atmosph ere, especially druin g 
the early phase of stellar evolution (jLammer et al.l[2C) 0Jl) , 
and energetic stellar cosmic rays could pose a serious 
risk to any form of surf ace life directly exposed to them 
(jGriefimeier et al.l[2005l ). 

The last step in order to estimate the magnetospheric 
properties and hence the level of magnetic protection is 
predicting the stellar wind properties for different stellar 
masses and as a function of planetary distance and time. 

There are two simple models used to describe the spa- 
tial structure and dynamics of the stellar wind: t he pure 
hydro dynamical model developed originally by iParked 
119581 that describes the wind as a non-magnetized, 
isothermal and axially symmetric flux of particles (here- 
after the Parker's model) and the more detailed albeit 
simpler niagneto- hydro dynami c model developed origi- 
nally by iWeber &: DavisI 119671 that takes into account 
the effects of stellar rotation and treats the wind as a 
magnetized plasma. 

It has been shown that Parker's model describes re- 
liably the properties of the stellar winds in the case of 
stars with periods of rotation of the same order of the 
present solar value, i.e. P ~ 30 days (|Preusse et al.l 
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l2005f ). However for rapidly rotating stars, i.e. young 
stars and/or active dM stars, the isothermal model un- 
deresti mates the ste llar wind properties almost by a fac- 
tor of 2 (jPreusse et al. 2005). For the purposes of scaling 
the properties of the planetary magnetospheres, (equa- 
tions ([3])-(l6|)), an underestimation of the stellar wind dy- 
namic pressure of that size, will give us values of the key 
magnetospheric properties that will be off by 10-40% of 
the values given by more detailed models. Magnetopause 
fields that have the largest uncertainties will be under- 
estimated by ^ 40%, while standoff distances and polar 
cap areas will be respectively under and overestimated 
by just ~ 10%. 

According to Parker's model the stellar wind average 
particle velocity v at distance d from the host star is 
obtai ned by solving the Parker's wind equation (jParkeJ 
[Ml): 



Using these formulae iGriefJmeier et al.ll2007al devised a 
clever way to estimate consistently T{t) and M^,(t) in the 
Parker's model and hence we are able predict the stellar 
wind properties as a function of d and t. For the sake 
of completeness we summarize here this procedure. For 
further details see section 2.4 in iGriefimeier et ani2007al 

For a stellar mass M^, and time t, the velocity of the 
stellar wind at d = 1 AU, uiau? is calculated using equa- 
tion (|37| . Replacing this velocity in the Parker's wind 
equation for d = 1 AU, we find numerically the temper- 
ature of the Corona T{t). This parameter is enough to 
provide us the whole velocity profile v(t, d, M^) at time 
t. To compute the number density we need the mass-loss 
rate for this particular star and at this time. Using the 
velocity and number density calculated from eqs. p7p 
and ([55)1 the mass-loss rate for the Sun Mq at time t 
and d — 1 AU can be obtained: 



u — log u = 4 log p-\ 3 

P 



(35) 



where u = v/vc and p = d/dc are the velocity and dis- 
tance normalized with respect to Vc — ^ksT /m and 
dc = GMi,m/{4kBT) which are respectively the local 
sound velocity and the critical distance where the stel- 
lar wind becomes subsonic. T is the temperature of the 
plasma which in the isothermal case is assumed constant 
at all distances and equal to the temperature of the stel- 
lar corona. T is the only free parameter controlling the 
velocity profile of the stellar wind. 

The number density n{d) is calculated from the veloc- 
ity using the continuity equation: 



n{d) 



M^ 



'i'Kd'^v{d)m 



(36) 



Here M* is the stellar mass-loss rate, which is a free 
parameter in the model. 

To calculate the evolution of the stellar wind we need a 
way to estimate the evolution of the coronal temperature 
T and the mass-loss rate M^,. 

Using observational estimates of the stellar mass-loss 
rate (Wood et al. 2002) and theoretical models for th e 
evolution of the stella r win d velocity ([N ewkirk 198(1) , 
IGriefimeier et al.l 120041 and iLammer etal.. .2004 devel- 
oped semiempirical formulae to calculate the evolution 
of the long-term averaged number density and velocity 
of the stellar wind for main sequence stars at a given 
reference distance (1 AU): 



viAu{t) = Wo ( 1 + - 



niAvit) = no ( 1 + - 



(37) 



(38) 



Here a^ 



, -0.4 3, a„ = 

(|GrieJ3meier et al.l [2009) 



-1.86 ± 0.6 and r = 25.6 Myr 

The parameters vq — 3971 

1.04 X 10^^ m~'^ are estimated from the 



km/s and no 

present long-term averages of the solar wind as measured 

at the distance of the Earth n(4.6Gyr, 1 AU, 1 M©) = 

6.59 X 106 m-3 and w(4.6Gyr, 1 AU, 1 Mq) = 425 km/s 

(|Schwennlll9"90l) . 



Meit) - 47r(l AU)2 m niAu(t)viAu(t) (39) 

Assuming that the mass-loss rate scales-up simply with 
the stellar surface area, i.e. M^{t) — MQ{t){Ri, / RqY , 
the value of M^, can finally be estimated. Using 
v{t^d^Mi,) and M* in the continuity equation (|36|) . the 
number density of the stellar wind n{t,d, Mi,) is finally 
obtained. 

The value of the stellar wind dynamic pressure 
Fdyn(i,d,M*) = mn{t,d,Mi,)v{t,d,A'Q^ inside the HZ 
of four different stars as computed using the procedure 
described before is plotted in figured 

It is important to stress here that for stellar ages 
i < 0.7 Gyr the semiempiric al formulae in eas. (l37ll 
and (|38l) are not longer reliable (jGrieCmeier et al.ll2007aD . 
These equations are based in the empirical relationship 
observed between the X-ray surface flux and the mass- 
loss rate M^, (jWood et all [2002 , 2005) which has been 
reliably obtained only for ages t > 0.7 Gyr. However 
iWood et al.l 120051 have shown that an extrapolation of 
the empirical relationship to earlier times overestimates 
the mass-loss rate by a factor of 10-100. At times t < 0.7 
Gyr and over a given magnetic activity threshold the 
stellar wind of mai n-sequence stars seems to be inhibited 
(jWood et al.ll2005l ) . Therefore the limit placed by obser- 
vations at t ~ 0.7 Gyr is not simply an observational con- 
straint but could mark the time where the early stellar 
wind also reaches a maximum (J.L. Linsky, Private Com- 
munication). This fact suggests that at early times the 
effect of the stellar wind on the planetary magnetosphere 
is much lower than normally assumed. Hereafter we will 
assume that intrinsic PMF are strong enough to protect 
the planet at least until the maximum of the stellar wind 
is reached at t sa 0.7 Gyr and focus on the stellar-wind 
and magnetosphere properties for times larger than this. 

6. RESULTS 

Using the results of our RTEM, the power-based scal- 
ing laws for dynamo properties, and the properties of 
the stellar insolation and stellar wind pressure, we have 
calculated the magnetosphere properties of earth-like 
planets and super-Earths in the HZ of different main- 
sequence stars. We have performed these calculations 
for hypothetical TPs in the mass-range 0.5-6 M^ and 
for the already discovered potentially habitable planets 
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GL 581d, GJ 667Cc and HD 40307g (see tabled]). The 
case of the Earth and an habitable Venus has also been 
studied for references purposes. 

To include the effect of rotation in the properties of 
the PMF we have assumed that planets in the HZ of 
late K and dM stars (M < O.TM©) are tidally locked 
at times t < 0.7 Gyr (n=2 in eq. ([55|. see figure [5]). 
Planets around G and early K stars (M > 0.7 Mq) will 
be assumed to have their primordial periods of rotation 
that we chose in the range 1 — 100 days as predicted by 
models of planetary formation (Miguel & Brunini 2010). 

Figures [7] and |8] show the evolution of magnetosphere 
properties for tidally locked and unlocked potentially 
habitable planets respectively. In all cases we have as- 
sumed that the planets are in the middle of the HZ of 
their host stars. 

Even at early times tidally locked planets of arbi- 
trary mass have a non-negligible magnetosphere radius 
Rs > 1.5 Rp. Previous estimates of the standoff dis- 
tances for tidally locked planets are much lower than the 
value s reported here. As an example iKhodachenko et all 
120071 place the standoff distances well below 2 Rp , even 
under mild stellar wind conditions (see figure 4 in their 
work) and independent of planetary mass and age. In 
contrast our model predicts standoff distances for tidally 
locked planets in the range of 2-6 Rp depending on plan- 
etary mass and stellar age. The differences between both 
predictions arise mainly from the underestimation of the 
dipole moment for slowly rotating planets found in these 
works. Thermal evolution and the dependency on plan- 
etary mass of the PMF properties are responsible for 
the rest of the discrepancies in previous estimates of the 
magnetosphere properties. 

Though tidally locked planets seem to have larger mag- 
netospheres than previously expected, they still have 
large polar caps, a feature that was previously over- 
looked. As a consequence of this fact well protected 
atmospheres, i.e. atmospheres that lie well inside of 
the magnetosphere cavity (hereafter magnetised planets) , 
could have more than 15% of their surface area exposed 
to open field lines where thermal and non-thermal pro- 
cesses could efficiently remove atmospheric gases. More- 
over, our model predicts that these planets would have 
multipolar PMF which contributes to an increase of the 
atmospheric a rea open to the interplan etary and magne- 
totail regions (Sisc oe fc Crookerlll976[ ). Then the expo- 
sition of magnetised planets to harmful external effects 
would be a complex function of the standoff distance and 
the polar cap area. 

Overall magnetic protection improves with time. As 
the star evolves the dynamic pressure of the stellar wind 
decreases more rapidly than the dipole moment (see fig- 
ures m and [6|) . As a consequence the standoff distance 
grows in time and the polar cap is shrunk. However with 
the reduction in time of the stellar wind pressure the 
magnetopause field is also reduced a fact that can affect 
the incoming flux of CR at late times. 

The sinuous shape of the contour lines in the middle 
and lower rows is a byproduct of the inner core solidifi- 
cation in planets with Mp < 2Mq. Critical boundaries 
between regions with very different behaviors in the mag- 
netosphere properties are observed at Mp ~ l.OM^ and 
Mp ^ I.8M0 in the middle and rightmost panels of the 
standoff radius and polar cap area contours. Planets to 



the right of these boundaries still have a completely liquid 
core and therefore produce weaker PMFs (lower standoff 
radius and larger polar cap areas). On the other hand, 
the inner core in planets to the left of the these bound- 
aries have already started to grow and therefore their 
PMFs are stronger. 

Unlocked planets (figure [S]) are better protected than 
slowly rotating tidally locked planets by developing ex- 
tended magnetospheres Rs ^ 4 i?p and lower polar cap 
areas Ape < 10%. It is interesting to notice that in 
both cases and at times i '^ 1 Gyr a smaller plane- 
tary mass implies a lower level of magnetic protection 
(lower standoff distances and larger polar caps) . This re- 
sult seems to contradict the idea that low-mass planets 
{Mp < 2) are better suited to de velop intense and pro- 
tective PMFs iG aidos et al. 2010: lTachinami et al.ll2011l : 
iZuluaga fc Cuarta£ll2012l ). To explain this contradiction 
one should take into account that magnetic protection as 
defined in this work depends on dipole moment instead 
of surface magnetic field strength. Since dipole moment 
scales-up as A^ ^ B^ipRp more massive planets will have 
a better chance to have large and protective dipole mo- 
ments. 

It is interesting to compare the predicted val- 
ues of the maximum dipole moment calculated 
here with the values roughly es ti mated in previous 
attempts (Gricfimcicr et al. '2005|; IKhodachenko et"al] 
2007: LoDCz-Moralcs ct al. 20111). On one hand 



Khodachenko et all 120071 estimate dipole moments for 
tidally locked planets in the range 0.022-0.15 M^. These 
values have been systematically used in the literature to 
study dif ferent aspects o f plane tary magnetic protection 
(see e.g. iLammer et al.l (|2010f) and references therein). 
For the same type of planets our model predicts max- 
imum dipole moments almost one order of magnitude 
larger (0.15-0.60 A^©) with the largest differences found 
for the most massive planets (M > 4Mq). These dif- 
ferences axis£jiran_rti£ja£tthatnone of the scaling-laws 
used bv IKhodachenko et al.l l2007 depend on the convec- 
tive power. In our results the dependency on power ex- 
plains the differences between massive planets and ligther 
planets especially at earl y times. On the other hand 
iLopez-Morales et al.ll201ll estimate magnetic dipolar mo- 
ments of tidally locked super-Earths in the range 0.1- 
1.0 Ai^. Thes e values are compatible wi th our results. 
In their, work Lo pez-Morales et al.l 12011! use the same 
power-based scaling laws we applied here but assume a 
rather simple interior model and a static thermal model 
where the convective power is set such that maximizes 
the efficiency with which the convective energy is con- 
verted into the magnetic field. 

A more detailed account of the evolution of magne- 
tosphere properties for the already discovered habitable 
planets is presented in figure [9] In all cases we have 
assumed that all planets have compositions similar to 
Earth (RTEM). Although almost all planets are tidally 
locked, we have also computed the magnetic properties 
for a primordial period of rotation P = 1 day. 

The case of the "hydrated" Venus is particularly in- 
teresting in order to analyse the rest of planets. The 
dynamo of the actual Venus probably shut down at 
t = 3 Gyr as a consequence of the drying of the man- 
tle (jChristensen et aTl |2009[ ) . A massive loss of water 



The influence of tfiermal evolution in the magnetic protection of terrestrial planets 



13 



Planet 


Mp(Me) 


RpiRfs) 


a(AU) 


Po (days) 


c 


S-typc 


Af*(Afo) 


agc(Gyr) 


tid. locked 


Rots. 


Earth 


1.0 


1.0 


1.0 


365.25 


0.016 


G2V 


1.0 


4.56 


No 


- 


Venus 


0.814 


0.949 


0.723 


224.7 


0.007 


G2V 


1.0 


4.56 


Probably 


- 


GJ 667Cc 


4.545 


1.5* 


0.123 


28.155 


< 0.27 


Ml. 25V 


0.37 


> 2.0 


Yes 


(1) 


GL 581d 


6.038 


1.6* 


0.22 


66.64 


0.25 


M3V 


0.31 


4.3-8.0 


Yes 


(2), (3) 


HD 40307g 


7.1 


1.7* 


0.6 


197.8 


0.29 


K2.5V 


0.77 


4.5 


No 


(4) 



TABLE 2 

Properties of the already discovered SEs inside the HZ of their host stars. For reference purposes the properties of 

Venus and the Earth are also included. Values of radii marked with an * are unknown and were estimated using the 

mass-radius relation for planets with th e sa me composition as th e earth, i.e. rj, = r (f)( mp/mff,)'^-^'^ (val06 ). references 

ARE: dl lBONFILS ET AL.|[20Tll. ('2') IUDRY ET AL.|[2007|. O") IMAYOR ET AL.ll2009l (4') ITUOMI ET AL.|[20Tg 



induced by a runaway greenhouse and insufficient early 
magnetic protection played a central role in the extinc- 
tion of the early Venusian PMF. The evolution of the 
PMF in the potentially habitable planets GL 581d, GJ 
667Cc and HD 40307g could have a similar fate. Their 
masses are much larger and therefore their atmospheres 
are protected by stronger gravitational fields. 

For planet HD 40307g our reference thermal evolution 
model predicts a late shut down of the dynamo i^yn ~ 4 
Gyr. According to our reference model the planet is 
presently devoid of a dynamo generated magnetic field. 
However, being around a K star {M^, ^ 0.7) the stel- 
lar wind and XUV radiation have probably decreased 
enough to not represent at present times a real threat 
for its atmosphere. 

GL 581d and GJ 667Cc are located in the HZ of dM 
stars where the stellar wind pressure and XUV radiation, 
even at times as late as 4 Gyr, are intense enough to erode 
their atmospheres or to make them lose their water con- 
tent. The RTEM predicts that for an Earth-composition 
GL 581d at present times had already lost its dynamo 
and has been exposed for almost 2.5 Gyr to the harmful 
effects of the stellar wind and CR. The planet however 
is the most massive of the three planets and probably 
has a thick atmosphere able to withstand the continuous 
agression of its host star. 

Given the estimated age of the GJ 667C system (i « 2 
Gyr), the RTEM predicts that the planet still has a dy- 
namo (red circle in figure [5]) . Magnetosphere properties 
are very close to that of our "hydrated" Venus, 4 Gyr 
ago. However its mass is lower than that of GL 581d 
and it is located at the inner edge of the HZ where the 
exposition to the XUV radiation from its host star (a 
young Ml star) could have been enough to induce mas- 
sive loss of atmospheric gases including water. We will 
come back on these issue in section WTU when we will show 
how to estimate the atmospheric mass-loss rate for this 
particular planet. 

6.1. Toward an estimation of the atmospheric mass-loss 

Combining the model of magnetosphere evolution de- 
veloped here with models of thermal and non-thermal 
atmospheric escape it would be possible to estimate the 
mass-loss rate from atmopsheres of magnetised and un- 
magnetised potentially habitable planets. This is a fun- 
damental goal to be pursued in the near future if we 
want to assess the actual habitability of present and fu- 
ture discovered terrestrial planets in the HZ of their host 
stars. The complex interaction between an inflated atmo- 
sphere and its protective magnetosphere and large uncer- 
tainties in the surface fluxes of atmospheric gasses that 
compensate the loss of volatiles induced by the action of 
the stellar wind, render this goal hard to achieve in the 



present. Despite these limitations we can still make order 
of magnitude estimations based on our own results and 
the mass -loss rate computed for exa mpl e in the recent 
work s bv lTian et a"Lll2008l lTianll2009l and lLammer et al.l 
[20T1 

Atmospheric thermal mass-loss induced by the ex- 
position to X-rays and EUV stellar radiation (XUV) 
have been est imated for the case of Earth-like N2 rich 
atmospheres (IWatson et al.l 119811: IKulikov et all 120061: 
iTian et all 120081) a.nd dry Venus-like C O2 rich atmo- 
spheres ( TianI l20"09t iLammer et al.l [20121) . One critical 
property of an inflated atmosphere is essential to evalu- 
ate the exposition of such atmospheres to further non- 
thermal processes: the radius of the exobase i?oxo- ^0x0 
is defined as the distance where the mean-free path of 
atmospheric particles could be comparable to the size of 
the planet. When the radius of the exobase is compara- 
ble or larger than the magnetic standoff distance Rs we 
will say that the planet is unmagnetised. Under these 
conditions the gases escaping from the exosphere will be 
picked-up by the stellar wind and lost to the interplan- 
etary space. On the other hand if the exobase is well 
inside the magnetosphere (which is the case of the Earth 
today) atmopsheric gasses escaping thermally from the 
exosphere could stay trapped by the magnetic field form- 
ing a plasmasphere. Planets under this condition will be 
magnetically protected and the mass-loss rate is expected 
to be much lower than for unmagnetised planets. 

Using the conservative estimation of the X and 
EUV luminosities of main-sequence stars given by 
iGarces et al.ll20lll we have estimated the XUV flux at the 
top of the atmospheres of GL 581d, GJ 667Cc and HD 
40307g during the first critical gigayear of planetary evo- 
lution. The planet that received the minimum amount of 
XUV radiation is HD 40307g with Fxvv = 35-10 PEV 
n PEV = 0.64 erg cni^^ s^^ is the Present Earth Value , 
[Judge et al.ll2003l : iGuinan et al.ll2009D . GL 581d was ex- 
posed in the first gigayear to a flux of Fxuv = 150 — 250 
PEV while GJ 667Cc received the maximum amount of 
XUV radiation among them, fxuv = 45 - 800 PEV. 

Using the recent results bv lTianI 120091 that computed 
the exosphere properties of massive super-Earths, i.e. 
Mp ^ 6M0, subject to different XUV fluxes, we can 
estimate the exosphere radius and mass-loss rate f or our 
three habitable super-Earths. Actually, since the iTianl 
(j2009( ) results are only available for planets with a min- 
imum mass of Mp = 6M0, a qualitative extrapolation 
of the results for 10, 7 and 6 M^ (see figure 4 and 6 in 
its paper), shows that exobase radius and mass-loss rates 
are larger for less massive planets. This is particularly 
useful at trying to apply the Tian's results to GJ 667Cc 
Mp ss 4.5Af0 and other less massive potentially habit- 
able planets. In these cases we will use the results by 
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iTianI (|2009[ ) to calculate a lower bound of the exobase 
radius and mass-loss rates. 

Using the estimated XUV flux for HD 40307g {Mp w 
7M9) and assuming an initial CO2 rich atmosphere, 
Tian's results predict that the exosphere of the planet 
and hence its mass-loss rate was low enough to avoid a 
significant early erosion of its atmosphere (see figures 4 
and 6 in his paper). This is true at least during the first 
1-2 Gyr during which our magnetic protection model pre- 
dict the planet was enshrouded by a protective magneto- 
sphere (see figure[9|). After dynamo shut down the atmo- 
sphere of HD 40307g has been exposed to the direct ac- 
tion of the steUar w ind. Assuming a stellar age of 4.5 Gyr 
(jTuomi et al.ll20l3) this eff'ect has been eroding the at- 
mosphere for 3-4 Gyrs. During this unmagnetised phase 
the atmospheric mass-loss rate can be simply estimated 
as M w amnvcs (IZendeias et al.|[20lol ) where a is the so- 
called entrainment efficiency and m, n and VeS are the 
mass, number density and effective velocity of the stellar 
wind as measured at planetary distance (see eq. ([2])). 
Using a entrainment efficiency a ~ 0.3 (which is appro- 
priate for example to describe the mass-loss of the Venus 
atmosphere) we obtain that the total mass-loss during 
the unmagnetised phase is less than 1% of a conser vative 
estim ate of the total volatile content of the planet (jTianl 
|2009( ). Although our model provides only upper limits 
to magnetic protection and the planet could have for ex- 
ample a lighter Nitrogen- rich atmos phere which is more 
prone to XUV induced mass-losses (iWatson et al.l 119811 : 
iKuhkov et all l2006t iTian et al.l |2008D, this preliminary 
estimation suggests that HD 40307g probably still pre- 
serve a dense enough atmosphere able to sustain surface 
liquid water and hence to be actually habitable. 

The case of GL 581d {Mp sa 6A/0) is quite different. 
Assuming that our estimations of t he XUV fiux are right, 
the exosphere radius predicted by iTianI 120091 should be 
close to the actual one. In this case at times i ^ 1 Gyr, 
^cxo = 1.8 — 2.3 i?p. However our reference magneto- 
sphere model predicts for this planet magnetic standoff 
distances Rs > 2.7 at all times. We conclude that using 
our estimations GL 581d could have been protected by its 
intrinsic magnetic field during the critical early phases of 
planetary evolution and probably has preserved the crit- 
ical volatiles in its atmosphere. Still the uncertainties 
in the exosphere model or in the atmopsheric composi- 
tion and of course in the magnetic model developed here 
should lead to a different conclusion and further theoret- 
ical and probably observational investigation is required. 

The most interesting case is that of GJ 667Cc {Mp « 
4.5M0). The minimum exosphere radius predicted for 
this planet at i ^ 1 Gyr lies between 3.0 — 4.5 Rp 
while according to our magnetosphere model, the mag- 
netic standoff distance is Rs < 3 Rp up to 2 Gyr. Since 
the exosphere radi us should ac tually be larger than that 
predicted with the iTianI (I2OO90 model and our magnetic 
model is actually optimistic, the chances that this planet 
was unprotected by its magnetic field in the critical first 
gigayear are high. But exposition does not necessarily 
mean a complete obliteration of the atmosphere (see for 
example the case of Venus) . To evaluate the level of ther- 
mal and non-thermal obliteration of the atmosphere we 
need to estimate the actual mass-loss rate. At the XUV 
fluxes estimated at the top of the atmosphere of this 



planet during the first gigayear, the minimum thermal 
mass-loss rate of Carbon atoms from a CO2 rich atmo- 
sphere will be larger than 2 — 4 x 10^° atoms cm~^ s~^ 
(see figure 6 in Tian 200^. We should recall that this is 
actually the value for a 6M(^ super-Earth. For the actual 
mass of the planet, 4.5M0, the mass- loss rate could be 
even larger. Moreover if as predicted here the exosphere 
is exposed directly to the stellar wind, non-thermal pro- 
cesses can contribute to a larger increase in the mass-loss 
from the planetary atmosphere. At the minimum mass- 
loss rate the exposed GJ 667Cc atmosphere could have 
lost more than ~ 10^^ atoms of Carbon in just ^ 100 
Myr and in the first gigayear the amount of carbon ther- 
mally lost to space could rise to ~ 10"*^^ atoms. If we 
scale-up linearly with planetary mass the total inventory 
of CO2 in the atmosphere, crust a nd mantle of Venus, 
which is 2 — 3 X 10''^ molecules fsee lTianl l2009 and refer- 
ences therein), a 4.5^^ planet will have a total budget 
of ~ 10^^ CO2 molecules. In summary at the minimum 
mass-loss rate and assuming a relatively rapid degassing 
of the planet, Gj 667C c could have lost its total inven- 
tory of Carbon to interplanetary space in the first couple 
of gigayears. Even assuming that large amounts of CO2 
are still trapped in the mantle and crust of the planet, 
its atmosphere should be being rapidly obliterated by 
the steUar wind. We speculate that GJ 667Cc is a sort 
of "Venus-like" planet. Regardless the fact the planet is 
well inside the radiative habitable zone it has lost its ca- 
pacity to support life via a massive stellar-wind induced 
loss of volatiles. 

7. DISCUSSION 

Applying simplified thermal evolution model and dy- 
namo scaling laws to planets whose bulk properties are 
barely known or even hypothetical it is challenging and 
probably raises more questions than it attempts to re- 
spond. Further observations of the potentially habit- 
able planets should be required to precise their physical 
properties and to reliably model its interior and thermal 
evolution. Moreover continued observational efforts to 
look for direct evidence or proxies of planetary magne- 
tospheres and any other signatures of magnetic protec- 
tion, though challenging, should also be attempted. Here 
we want to discuss the assumptions on which our global 
model relies, its uncertainties as measured by the sensi- 
tivity of the model to changes in key physical parameters 
and the missing pieces of information and present obser- 
vational limitations to confirm or improve this and other 
models of planetary magnetic protection. 

7.1. Model assumptions 

The strength of a physical model depends on the hy- 
pothesis and assumptions on which it relies. Apart from 
numerous albeit very common assumptions, the magnetic 
protection model presented here depends on three major 
assumptions we discuss in the following paragraphs. 

We have assumed that terrestrial planets always de- 
velop an initially metallic liquid core, irrespective of 
their composition and early formation history. This 
is not necessarily true. The formation of a metallic 
liquid core would depend on very complex processes 
and other barely known physical factors. It has been 
shown for example that under extreme water oxidation 
of iron the formation of a metallic core will be avoided 
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(|Elkins-Tanton fc Seageil2008f) . In this case silicate core- 
less planets will be formed. On the other hand even if a 
planet is well differentiated the core could be solid from 
the beginning (see e.g. VAL06). However, it should be 
emphasized here that our model provides only the best- 
case scenario of magnetic protection. Therefore, if under 
the assumption of having a liquid metallic core, a planet 
is found to be lacking enough magnetic protection, the 
case when the planet is not well-differentiated or when 
it never develops a liquid core will be even worse. In 
these cases the conclusions drawn from our model will 
be unchanged. 

The calculation of key magnetosphere properties re- 
lies on very simplistic assumptions about the complex 
physics behind the interaction between planetary and in- 
terplanetary magnetic fields and stellar wind. In partic- 
ular, standoff distances calculated with eq. ([5]) assume 
a negligible plasma pressure inside the magnetosphere. 
This condition could be violated in planets with very in- 
flated atmospheres and/or at close-in orbits. Under these 
extreme conditions the magnetic standoff-distance given 
by eq. ^ could be a poor underestimation of the actual 
magnetopause distance. However, the weak dependence 
of standoff distances and polar cap areas of the stellar 
wind dynamic pressure, offers some idea as to the ac- 
tual role that magnetospheric plasma pressure may pay 
in determining the size of the magnetosphere. Adding a 
plasma pressure term Ppj to the magnetic pressure P^p 
in the left-hand side of eq. ([3]), is equivalent to sub- 
stracting it from the stellar wind dynamic pressure Pgw 
in the right-hand side. An effective stellar wind pressure 
Pg'^ — Psw(l — Pp\/Psw) would replace the stellar wind 
term in the standoff distance definition (eq. (O). As a re- 
sult a plasma-pressure correcting factor (1— Ppi/Psw)^^^^ 
will modify our estimated purely-magnetic standoff dis- 
tance. Even in a case where the plasma was able to exert 
a pressure 50% of that of the stellar wind, the standoff 
distance will be increased by only a 10%. On the other 
hand in order to have a standoff distance one order of 
magnitude larger than that estimated using eq. ([5|) the 
plasma pressure inside the magnetosphere should amount 
for 99.999% of the total pressure. This is precisely what 
an unmagnetised planet would look like. In summary in- 
cluding more realistic condition into the definition of the 
magnetosphere boundary will not modify too much our 
results. 

A final but not less important assumption in our model 
is that of quiet stellar wind conditions. We have only 
taken into account average or quiet stellar wind condi- 
tions. We have completely neglected the effects of large 
but transient conditions such as those produced during 
coronal mass ejections (CME). To model the effect that 
a steady fiux/influx of CME plasma could have in close- 
in planets, we can modify the stellar wind pressure by 
maintaining the nominal velocity of the plasma but in- 
creasing the number density of wind particles by a factor 

of twc0. Taking into account that Rs ~ J^sw ^^ found 
that under the harsher conditions the magnetosphere ra- 
dius and polar cap areas will be modified only by 10-30% 
in respect to the nominal or quiet values computed here. 

® Under typical conditions of solar CME, the velocity of the wind 
is not modified too much but the plasma densities is increased up 
to 5-6 times over the average particle density 



This simplified estimation shows that our results seem to 
be robust in relation to uncertainties in the stellar wind 
pressure. However given the complexity of the interac- 
tion between the magnetosphere and the stellar wind un- 
der active phases, a further examination of this case is 
required and is left to future research. 

7.2. Sensitivity analysis 

In order to study the effect that uncertainties in sev- 
eral critical thermal evolution parameters have in the 
prediction of the overall magnetic protection of poten- 
tially habitable TPs, we performed a sensitivity analysis 
of our model. For this purpose we independently varied 
the value of 6 carefully chosen parameters of the model 
(see below) and compared the predicted dipole moment, 
the time of inner core formation and the dynamo lifetime 
with the same values obtained using the RTEM. 

We performed these comparisons for planets with five 
different masses: 0.7, 1.0, 3.5, 4.5 and 6.0 M^ (see fig- 
ure nop . These masses correspond approximately with 
those of the already discovered habitable planets (see ta- 
ble [2]) including a hydrated Venus and present Earth. In 
all cases we assume for simplicity a primordial period of 
rotation of P = 1 day. 

Since the dipole moment is an evolving quantity we 
have plotted in figure [TU] the average value of this quan- 
tity as calculated in the interval 0.7-2.0 Gyr. For times 
earlier than 0.7 Gyr the stellar wind pressure is uncer- 
tain and the magnetic protection cannot be estimated (as 
discussed in section 15.21 observations suggest lower stel- 
lar wind pressures at times earlier than this). For times 
larger than 2.0 Gyr the fiux of XUV radiation and the 
stellar wind pressure has decreased below the initial high 
levels. Although an average of the dipole moment is not 
phenomenologically relevant, it could be used as a proxy 
of the overall magnetic shielding of the planet during the 
harsh early phases of stellar and planetary evolution. 

After studying the full set of physical parameters in- 
volved in our interior structure and thermal evolution 
models (see table [T]) we identified 6 critical parameters 
whose values could have noticeable effects on the results 
or are subject to large uncertainties. We performed an 
analysis of the sensitivity that the model have to the 
variation of the following physical parameters: 

1. The core mass fraction, CMF. This is the 
fraction of the planetary mass represented by the 
metallic core. This parameter is determined by the 
Fe/Si ratio of the planet that it is fixed at planetary 
formation or could be altered by exogenous pro- 
cesses (e.g. late large planetary impacts). Our ref- 
erence model uses the Earth's value CMF — 0.325, 
i.e. assumes that all planets are dominated by a 
sillicate-rich mantle. As a comparison Mars has a 
CMF = 0.23 and the value for Mercury is CMF = 
0.65 (it should be recalled that Mercury could have 
lost a significant fraction of its mantle sillicates in- 
creasing the total iron fraction, probably after an 
early large impact). The CMF determines the size 
of the core and hence the thermal properties of the 
convective shell where the magnetic field is gener- 
ated. For our sensitivity analysis we have taken 
two extreme values of this parameter, CMF = 0.23 
(a mars-like core) and CMF = 0.43 (an iron-rich 
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core). Planets with larger cores have low pressure 
olivine mantles and our rheological model becomes 
unreliable. 

2. The initial temperature contrast at the 
CMB, ATcMB = eadbATadb (see eq. p|)). This 
is one of the most uncertain properties in thermal 
evolution models. The initial core temperatures 
would be determined by random processes involved 
in the assembly and differentiation of the planet. It 
could vary widely from planet to planet. In order 
to fit the thermal history of the Earth (time of in- 
ner core formation, present size of the inner core 
and magnetic field strength) we set Cadb = 0.7 and 
applied the same value to all planetary masses. In 
our sensitivty analysis we varied this parameter be- 
tween two extreme values of 0.6 and 0.8. Though 
we are not sure that this interval is representative 
of planets with very different masses and formation 
histories, our analysis provides at least the magni- 
tude and sign of the effect that this parameter has 
in the dynamo properties predicted by our thermal 
evolution model. 

3. High pressure viscosity rate coefficient, b (see 

eq. (|26l) ). Rheological properties of sillicates inside 
the mantle are among the most uncertain aspects 
of thermal evolution models. They critically de- 
termine, among other key quantities, the amount 
of heat that the core and mantle could transport 
through their respective boundary layers (see eqs. 
dni), dSni) and (121])). We found that the viscosity 
of the lower mantle (perovskite) is the most impor- 
tant source of uncertainties in our thermal evolu- 
tion model. The formula used to compute viscos- 
ity at that layer (see eq. ([27| ) strongly depends on 
temperature and pressure and the parameter con- 
trolling this dependence is the "rate coefficient" b. 
In the RTEM we used the value b = 12.3301 to re- 
construct the thermal properties of the Earth. In 
this value all the figures are significative reflecting 
the strong sensitivity of the model to this parame- 
ter. To study the impact of b in the model results, 
we varied it in the interval 10 — 14. 

4. Iron thermal conductivity, kc- This param- 
eter controls the amount of heat coming out 
from the core. In the RTEM we used a value 
A;c=40 Wm^-'^K^^ that fits the thermal evolution 
history and present magnetic field of the Earth 
(see table [1]). Although recent first-principles 
analysis suggest that values as large as 150-250 
Wm~^ K~^ could be comm on at Earth's core con- 
ditions (jPozzo et al.ll2012l ) we conform here to the 
standard values of this parameter. Further inves- 
tig ations to explore values as large as that found 
by iPozzo et al.ll20l3 should be attempted. For our 
sensitivity analysis we varied kc between 35 and 70 
Wm~^K~^, two values which are inside the typical 
uncertainty assumed for this property. 

5. Griineisen parameter for iron, 70c. This is one 
of the most critical parameters of the equation of 
state specially at core conditions. It affects strongly 
the mechanical structure, temperature profile and 



phase of iron in the metallic core (for a detailed 
discussion on the sensitivity of interior structure 
models to this parameter see e.g. VAL06). In the 
RTEM we used the reference value 7oc=1.36 that 
fits the thermal evolution history and present mag- 
netic field of the Earth (see table [Ij. Assuming dif- 
ferent kind of core alloys a relatively large range of 
values of this parameter has been used in literature 
(see VAL06 and references there in). Griineisen 
parameter values have been found in the range of 
1.36-2.338. Since our RTEM value is at the lower 
end of this range for our sensitivity analysis we have 
recalculated the model for a larger value of 2.06. 

Other uncertain parameters such as the critical 
Rayleigh number at the CMB, Ruc, that is also var- 
ied to study the sensitivity of thermal evolution models 
(i Gaidos. E]l2011h . were also studied. We did not found 
significant sensitivity of our model to variation of those 
parameters. 

We depict in figure [10] the relative variation in the 
aforementioned magnetosphere and dynamo properties 
when each of the previously described parameters were 
varied independently. 

We have found that planetary composition (CMF), 
mantle viscosity are responsible for the largest uncertain- 
ties in the predicted magnetic properties of the planet. 
Planets with small metallic cores have on average low 
magnetic dipole moments (squares in the first column of 
the upper panel). This is mainly due to a geometrical 
effect. The total heat produced by the core and hence 
the magnetic field strength at core surface is of the same 
order for Fe-poor and Fe-rich planets. However a small 
core means also a lower magnetic dipole moment, i.e. 
Ai ^ Re- Planets with lower content of iron also have 
small and hot cores and therefore the solid inner core for- 
mation and the shutting down of the dynamo are slightly 
delayed (squares in the middle and lowest panel of figure 

do]). 

Viscosity dependence on pressure and temperature, as 
quantified by the parameter 6, has the opposite effect 
on planetary magnetic properties than CMF at least for 
earth-like planets. A low viscosity lower mantle will 
favour the extraction of heat from the metallic core in- 
creasing the available convective energy for dynamo ac- 
tion. On the other hand a viscuous lower mantle will 
delay the formation of a solid inner core and extend the 
lifetime of the dynamo (middle and lowest panel in figure 

The effect of the Griineisen parameter at core condi- 
tions are negligible, at least on what respect to the mag- 
netic field strength and lifetime (upper and lower panels) 
which are the most critical properties affecting planetary 
magnetic protection. Only the time of inner-core forma- 
tion is strongly affected by changes in this parameter. 
In planets smaller than Earth, inner-core solidification 
can be delayed up to three times the reference value. On 
the other with a larger Griineisen parameter hand, earth- 
like planets could get a solid inner-core very early in their 
thermal histories even almost since the beginning. This 
is the result of the interplay between the resulting evo- 
lution of the thermal profile and the solidus. 

The effect of the initial temperature contrast across 
the CMB, quantified by the parameter eadb, goes in the 
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same direction as viscosity. The reasons for this behavior 
are however far more complex. A larger initial tempera- 
ture contrast across the CMB also implies a larger initial 
temperature at the core center. Although a hotter core 
also produce a larger amount of available convective en- 
ergy, the time required for iron to reach the solidification 
temperature is also larger. The dynamo of planets with 
Mp < 2Mq and hot cores (large eadb) is weaker than that 
of more massive planets during the critical first couple of 
gigayears where the average is calculated. Planets with 
a colder core develops a solid inner core almost from the 
beginning and the release of latent and gravitational en- 
ergy feeds a stronger dynamo. 

The case of more massive planets, Mp > 2M^, where 
the condition for an inner core formation is never reached 
during the dynamo lifetime, is different. In this case 
planets with hot cores (large lower mantle viscosities 
or high temperature contrasts along the CMB) produce 
large amounts of available convective energy. A larger 
convective power will produce a larger value of the local 
Rossby number. Thus massive planets with hot cores 
also have multipolar dynamos and hence lower dipole 
magnetic moments and a reduced magnetic protection. 

Thermal conductivity affect less the results of the ther- 
mal evolution model. Besides the case of massive planets 
where differences in the order of 10-30% in the magnetic 
properties are observed when kc varies between its ex- 
tremes, the magnetic properties calculated with our ref- 
erence thermal evolution model seem very robust against 
variations in these two properties. However, it should be 
mentioned that this result applies only when a standard 
value of kc is assume d. Further investig ations to explore 
the recent findings (jPozzo et al.l 120121 ) Concerning the 
possibility that kc could be larger by factors as large as 
2-3 should be attempted. 

In summary, despite the existence of a natural sensitiv- 
ity of our simplified thermal evolution model to uncer- 
tainties on their free parameters, the results presented 
in this paper seem to be correct at least in the order of 
magnitude. Moreover since the standoff distance and po- 
lar cap areas, which are the actual proxies to magnetic 
protection, goes as M^^^, a one order of magnitude es- 
timation of Ai will give us an estimation of the level of 
magnetic protection off by a factor around 2. 

To clarify this point let's consider the case of GL 581d. 
If we assume for example that its iron content is much 
less than that of the Earth (as was assumed in the RTEM 
model) but the rest of critical thermal properties are es- 
sentially the same, the average standoff distance (polar 
cap area) at the critical first gigayear will be off by only 
^ 30% with respect to the prediction depicted in figure[2l 
More interesting is to notice that probably the uncertain- 
ties due to the unknown period of rotation (shaded area 
in figure [9]) seem to be much larger than those coming 
from the uncertainties in the thermal evolution model. 

7.3. Observational support 

Validating or improving thermal evolution models and 
dynamo scaling laws for the case of super-Earths nowa- 
days represents a huge observational challenge. The 
available sensitivity of our best instruments in the Earth 
and in space are rather insufficent. New and/or im- 
proved instruments and observational techniques will 
be required to detect, catalogue and compare the ther- 



mal and magnetic properties of low mass planets in the 
medium to far future. However the importance that the 
detection and characterization of the magnetic proper- 
ties of future discovered potentially habitable planets in 
order to assess their true habitability clearly justify the 
effort. 

The first goal seems to be the direct or indirect de- 
tection of super-Earth magnetospheres. Four methods, 
some of them already used in our own solar system, could 
be devised to achieve this goal: 1) the detection of radio 
waves coming from synchrotron and cyclotron radiation 
produced by plasma trapped in the magnetosphere, 2) 
the detection of a bow shock or a tail of ions produced 
by the interaction of the planetary atmosphere and mag- 
netosphere with the stellar wind or the interplanetary 
plasma, 3) the detection of planetary auroras and 4) spec- 
troscopic observations of a non-equilibrium atmospheric 
chemistry induced by a high flux of CR (this is actually 
a negative detection of a magnetosphere) . 

The first (radioemission) and second methods (bow 
shock or tail) have already b een studied in detail 
jBastian et al. 2000"; 'Farrell et al."2004' 'Griefimeier et al.| 
|2007b; Lazio et al. 2009; Vidotto et al. 2010, 2011.). Its 
reliability, at least for the case of planets with intense 
magnetic fields or placed very close to their host star, 
has been already tested. 

If synchrotron or cyclotron radiation coming from the 
magnetosphere of terrestrial planets could be detected, 
the power and spectra of the radiation could be used 
to measure the magnetic field strength. However, even 
with the most sensitive instruments, e.g. the Low Fre- 
quency Array, LOFAR or the Long Wavelength Array 
LWA, the expected power and spectra are several orders 
of magnitude below the threshold of detection. Pow- 
ers as large as 10^ - 10^ times the Jupiter radioemission 
and frequencies in the range of tens of MHz are required 
for the present detection of synchrotr on radio emission 
in plan etary magnetospheres (see e.g. iGriefimeier et al.1 
l2007bD . The magnetic field intensities and expected fre- 
quencies produced in super-Earths magnetospheres are 
several orders of magnitude lower than these thresholds 
and probably are far from being detected in the near 
future. 

It has been shown recently that measurements of the 
asymmetry in the ingress and egress of transiting planets 
can be used to detect the presence of a bow shock or a tail 
of plasma around the planet. Vidoto (2010, 2011) used 
this phenomenon to constrain the magnetic properties 
of Wasp-12b. The formation of a detectable bow-shock 
depends, among other factors, on the relative velocity be- 
tween the planet and the shocked plasma. Close-in plan- 
ets with strong enough magnetic fields (this is precisely 
the case of Wasp- 12b) can easily develop UV-opaque bow 
shocks and allow reliable detection. However low-mass 
planets with relatively weak magnetic fields such as those 
predicted with our models, hardly produce a detectable 
bow-shock. It has been estimated that magnetopause 
fields in the range of several Gauss should be required to 
have a detectable signal of a bow shock (A. A. Vidotto, 
Private Communication). Our habitable super-Earths 
have magnetopause fields in the order of a few micro 
Gauss (see figure [7]). The case of an ion tail coming from 
a weakly magnetised planet has received less attention 
and probably could offer better chances for a future indi- 
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rect detection of the magnetic environment of low mass 
planets. 

Finally the detection of far UV (FUV) or X-ray emis- 
sion from planetary auroras can also be used as a tool to 
study directly and indirectly planetary magnetospheres. 
Planetary auroras with intensities as high as 10^ to 10"^ 
times larger than that of the Earth are expected in close- 
in giant p lanets subject to t he effects of CMEs from its 
host star (|Cohen et al.ll2Qll"l ). If we estim ate that a typi- 
cal Earth Aurora has an intensity of 1 kR (iNeudegg et alJ 
I200H) (Being 1 R - 10"" photons m-^s-igrad-^) and 
assuming that 10% of a close-in Jupiter-like planet is 
covered by auroras producing FUV photons around 130 
nm, the total emitted power from these planets will be 
r^ 10^^ W. If we assume that this is the present thresh- 
old for exoplanetary aurora detection, even under strong 
stellar wind conditions and distances typical of the hab- 
itable zone, the total FUV power produced by auroras 
in the polar cap of earth-like planets could be only 10^ 
W which is 8 orders of magnitude less than the present 
detection threshold. Not to mention that the FUV ra- 
diation should be detected against an intense UV back- 
ground coming from a probably young and active low 
mass star. If we can find ways to overpass these diffi- 
culties, the observation of the FUV and X-ray emmision 
and its variability from auroras in potentially habitable 
super-Earths could be used as powerful probes of the 
magnetic environment around the planet. 

8. SUMMARY AND CONCLUSIONS 

We studied here the influence that the thermal evo- 
lution of potentially habitable terrestrial planets has in 
the protection that an evolving planetary magnetopshere 
could provide against the atmospheric erosion caused by 
the stellar wind. 

We developed a simple parametrized thermal evolu- 
tion model able to reproduce the global thermal history 
and magnetic properties of the past and present Earth. 
We applied this model to predict the thermal histories of 
planets with masses in the range of 0.5 to 6.0 AIq and 
with chemical compositions similar to Earth. Using these 
results and applying up-to-date dynamo scaling laws we 
predicted the magnetic properties of terrestrial planets in 
the habitable zone as a function of time, planetary mass 
and rotation rate. A simple model of the evolution and 
interaction of the stellar wind with the planetary mag- 
netic field, that has been adapted from previous works, 
allowed us to compute the global properties of the mag- 
netosphere in order to assess the level of magnetic protec- 
tion that potentially habitable Earth-like planets could 
actually have. 

We applied our model to the case of already known 
potentially habitable terrestrial planets (GL 581d, HD 
40307g and GJ 667Cc), to the Earth itself and to the case 
of a hydrated Venus. In the case of the Earth our model 
reproduces fairly well the early and present thermal and 
magnetic properties of our planet. In the case of the hy- 
drated Venus, the model predicts low values the standoff 
distance and large polar cap areas in the frist critical 
gigayear of planetary and thermal evolution, which are 
compatible with the idea that the planet lacked a strong 
enough magnetic protection able to avoid a massive loss 
of water and volatiles that finally lead to the shut down 
of its dynamo ~ 3 Gyr ago. 



Compelling results were found in the case of the three 
already discovered extra-solar-system potentially habit- 
able planets. Assuming an earth-like composition and 
thermal evolution parameters similar to those used in 
the case of the Earth (reference thermal evolution model, 
RTEM), our model predicted that the dynamo of GL 
581d and HD 40307g have been already shut down. A 
younger GJ 667Cc seems to still have an active dynamo. 

A non-trivial dependence of the magnetic properties on 
planetary age, planetary mass and period of rotation has 
been found in general for terretrial planets inside the HZ 
of their host stars. Thermal evolution is responsible for 
the non-trivial relationship among all these properties. 
Contrary to what was found in previous work, tidally 
locked planets could develop relatively intense magnetic 
fields and extended magnetospheres. However they also 
have extended polar caps and probably multipolar mag- 
netic fields where field lines open to the interplanetary 
space and magnetotail regions probably increasing the 
non-thermal mass-losses. 

Using recent results for the relationship between the 
exposition to XUV radiation, the exobase radius and 
mass-loss rate from massive super-Earths, we estimated 
the level of exposure and mass- losses for the three already 
discovered potentially habitable super-Earths. With the 
available information not too much could be said about 
the magnetic protection of HD 40307g. Further theo- 
retical investigations are required to evaluate this case. 
Our model predicts a large enough magnetosphere able to 
protect GL 581d against the erosive action of the stellar 
wind during the first critical phases of planetary evolu- 
tion. However since our model is still optimistic further 
theoretical and probably observational analyses should 
be performed to establish on a more solid basis the mag- 
netic protection of this planet. Our upper-limit to the 
standoff-distance and the most optimistic estimation of 
exobase radius and mass-loss rate from the atmosphere 
of GJ 667Cc, point-out the fact that this planet has al- 
ready lost a large fraction of its inventory of volatiles. 
All the evidence compiled in this work make GJ 667Cc a 
sort of "Venus analogue". Although further theoretical 
analyses are required our best guess is that despite the 
fact that it is inside the radiative HZ of its host star the 
planet is presently uninhabitble. 

Last but not least we tested the robustness of our con- 
clussions by changing several of the most sensitive input 
parameters of our thermal evolution model. We found 
that even under the present uncertainties the predicted 
properties of planetary magnetopsheres are rather ro- 
bust. We calculated that introducing large variations in 
the composition of the planets and the rheological and 
thermal properties of their interiors with respect to the 
reference thermal evolution model, the critical magnetic 
properties, such as the standoff radius and the area of 
the polar cap, change only by a factor of two. Results 
are also robust against uncertainties in the stellar wind 
properties that could be very important in the case of 
close-in habitable planets around active and young dM 
stars. 

The problem of evaluating the magnetic protection of 
potentially habitable planets is far from being settled. 
Other sources of intrinsic magnetic fields, thermal evolu- 
tion and interior structure of planets with "exotic" com- 
positions, improved theoretical models and new experi- 
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mental evidence of the behavior of iron at high pressures 
and temperatures, improved and validated models of the 
evolution and spatial structure of stellar winds and of 
course more and better observational data coming from 
the already discovered habitable super-Earths and future 
discovered potentially habitable exoplanets, will allow us 
to assess the actual magnetic protection of potentially 
habitable environments. 
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Fig. 1. — Schematic representation of the planetary interior. In the schematic sUce we depict the main quantities used here to describe 
the thermal evolution of the planets. The temperature profile depicted below the slice does not use real data. Distances and sizes are not 
represented with the right scale. 
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Fig. 6. — Evolution of the stellar wind dynamic pressure at the center of the HZ for a selected set of stellar masses. The reference average 
solar wind pressure is Psw© = 1-86 nPa. Dashed curves indicate the value of the stellar wind pressure at the inner and outer edges of the 
HZ around stars with 0.2 Mq and 1.0 Mq respectively. The HZ limits where the pressure were calculated are assumed static and equal to 
those at T = 3 Gyr. Stellar wind pressure at i < 0.7 Gyr computed with the semiempirical model used in this work is too uncertain and 
were not plotted. 
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Fig. 7. — Evolution of the magnetopause field (upper row), standoff distance (middle row) and polar cap area (lower row) of tidally 
locked (slow rotating) planets around late dK and dM stars. The rotation of each planet is assumed equal to the orbital period at the 
middle of the HZ (see values in the rightmost vertical axes) . The value of the magnetosphere properties rerturned by the contour lines in 
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Fig. 9. — Evolution of magnetosphere properties for the already discovered habitable SEs, the Earth and an "hydrated" version of Venus 
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rotation period of theplanet. Filled circles are the predicted present day magnetosphere properties computed according to the properties 
summarized in table [2] 
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Fig. 10. — Sensitivity analysis of our reference thermal evolution model (RTEM). Squares and diamonds indicate the relative value of 
three critical magnetic properties, {Mdip ) (the average of the dipole moment between 0.7 and 2 Gyr), tic (time of inner core formation) 
and tdyn (dynamo lifetime), as calculated by the thermal evolution model. For the analysis 5 different key thermal evolution parameters 
were independently changed with respect to the reference value in the RTEM: the core mass fraction (CMF), the thermal conductivity of 
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